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1.1 Abstract 



Two models involving particles moving by "hopping" in disordered media 
are investigated: 

A model glass-forming liquid is investigated by molecular dynamics un- 
der (pseudo-) equilibrium conditions. "Standard" results such as mean 
square displacements, intermediate scattering functions, etc. are reported. 
At low temperatures hopping is present in the system as indicated by a 
secondary peak in the distribution of particle displacements during a time 
interval t. The dynamics of the model is analyzed in terms of its potential 
energy landscape (potential energy as function of the 3N particle coordi- 
nates), and we present direct numerical evidence for a 30 years old picture 
of the dynamics at sufficiently low temperatures. Transitions between local 
potential energy minima in configuration space are found to involve particles 
moving in a cooperative string-like manner. 

In the symmetric hopping model particles are moving on a lattice by 
doing thermally activated hopping over energy barriers connecting nearest 
neighbor sites. This model is analyzed in the extreme disorder limit (i.e. 
low temperatures) using the Velocity Auto Correlation (VAC) method. The 
VAC method is developed in this thesis and has the advantage over previous 
methods, that it can calculate a diffusive regime in finite samples using 
periodic boundary conditions. Numerical results using the VAC method 
are compared to three analytical approximations, including the Diffusion 
Cluster Approximation (DC A), which is found to give excellent agrement 
with the numerical results. 
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Chapter 2 

Introduction 



This Ph.D. thesis deals with two models characterized by particles mov- 
ing by "hopping" in disordered media. The first model is a simple model 
glass-former, investigated here by computer simulations using molecular dy- 
namics [|J|. Glass forming liquids are liquids that upon cooling falls out of 
equilibrium and form a glass, see figure 2A. For general reviews on glass- 
forming liquids see H-H. For reviews on computer simulations of glass- 
forming liquids see (ft |ll| . 




Liquid 



Figure 2.1: Schematic overview of the behavior of glass-forming liquids. 
Upon cooling crystallization is avoided and the liquid becomes "super- 
cooled", which is a (pseudo-) equilibrium state (the crystal being the real 
equilibrium state). Upon further cooling the liquid falls out of equilibrium 
and undergoes a glass-transition at the temperature T g . T g depends on 
cooling rate; glass II is obtained by a lower cooling rate than glass I. 
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CHAPTER 2. INTRODUCTION 



The normal liquid behavior is to crystallize upon cooling below the melt- 
ing temperature, Tl- In a large class of liquids the crystallization can be 
avoided (e.g. by fast cooling) and the liquid becomes "supercooled". A 
supercooled liquid is in a pseudo-equilibrium state; The crystal is the real 
equilibrium state, but the supercooled state is typically stable on very long 
time scales. The supercooled state is characterized by relaxation times in- 
creasing strongly with decreasing temperature. If cooled with a constant 
cooling rate this means that the liquid at some temperature falls out of 
equilibrium and undergoes a glass-transition and forms a glass, which is a 
disordered solid. The glass-transition temperature, T g , depends on the cool- 
ing rate; a lower cooling rate results in a lower glass transition temperature. 
By convention the "laboratory" glass transition is taken to be where the 
relaxation time r is on the order of 100 sec (i.e. the "typical" time scale in 
the laboratory )[]. 

Binary systems with particles differing about 20% in diameter have been 
found to be good candidates for simple models of glass-forming liquids in 
computer simulations (8|,p, 12-ll|. At sufficiently low temperatures particles 



in these systems are found to move by "hopping" , i.e. they are localized for 
a period of time, and then move more or less directly to another position, 
where they again become localized [||,[9, 12 |16| . Since hopping of some form is 



expected to play a increasingly dominant role as temperature is lowered fl9| , 
pC|] , systems exhibiting hopping at temperatures which can be reached under 
equilibrium conditions in computer simulations are of special interest. In this 
thesis we investigate the dynamics of a binary Lennard- Jones liquid, which 
has earlier been shown to exhibit hopping [|n|. This is done under (pseudo-) 
equilibrium conditions, i.e. above the glass transition temperature. 

The model glass-former investigated in this thesis is not a model of 
any particular liquid existing in the laboratory (or the real world for that 
matter). Like most models glass formers studied in computer simulations, 
it should be thought of as a (very) simple model exhibiting some of the 
complex behavior found in real glass forming liquids. The reason why it 
is interesting to investigate such a simple model is two-fold: I) The glass- 
transition is found to occur in liquids that are chemically very different, but 
the related phenomenology is found to be strikingly similar. Since liquids 
with different chemical details behave in a similar way, it makes sense to 
study simple models ignoring the chemical details, in an attempt to un- 



1 The relaxation time is related to the viscosity r) by r\ — G x t, where Goo is the 
instantaneous shear modulus, r ~ 100s typically corresponds to r\ ~ 10 13 poise, which is 
sometimes used as the definition of the laboratory glass transition. 
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derstand the fundamental mechanisms. II) Most theories for glass-forming 
liquids also treat these in a highly simplified manner, ignoring the chemical 
details (for the same reasons given above). Consequently computer simula- 
tions makes it possible to test these theories on their own terms, i.e. without 
any additional approximations. One particular example of this is the Mode 
Coupling Theory (MCT) |21, 22|. When doing computer simulations of glass 
forming liquids it is rather natural to compare the results with the predic- 
tions of the MCT since these are plenty and detailed, and MCT deals with 
the temperature range accessible (at the present) to computer simulations 
under equilibrium conditions. It is not a goal in it self to test MCT in the 
present thesis, but some of the results are compared to the predictions of 
the MCT. 

The second model investigated in the present thesis is a so-called "hopping- 
model", where the hopping behavior of particles is "built in"; In the sym- 
metric hopping model |23-2^| particles are moving on a lattice by doing 
thermally activated hopping over random energy barriers. This model is 
found to have interesting universal features in the extreme disorder limit, 
i.e. when the temperature goes to zero. The symmetric hopping model 
has earlier been treated as a model for frequency dependent conduction in 
glasses p^,|2a|, with particles representing non- interacting charge carriers 
(ions or electrons). In that context the universality manifests itself as fol- 
lows: At low temperatures the shape of the conduction vs. frequency curve 
becomes independent of temperature, which is also what is found in experi- 
ments (see eg. |29|-31|). In the present thesis the symmetric hopping model is 
treated in the slightly more general context of diffusion in disordered media. 
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CHAPTER 2. INTRODUCTION 



2.1 Outline 

This thesis consists of 3 main chapters: 

In chapter 3 we report results from the simulations of the model glass- 
former mentioned above, with emphasis on the dynamics. In this chapter 
"standard" results are reported, i.e. the measures commonly used to de- 
scribe (glass-forming) liquids. All of the features found for this particular 
model has been reported earlier for other models, and this chapter thus con- 
stitutes a "review" of the behavior of model glass-formers, illustrating the 
similarities in the behavior of these. 

In chapter 4 the dynamics of the model glass-former described in chapter 
3 is investigated in term of its "potential energy landscape" , which is simply 
the potential energy as a function of the 3iV particle coordinates. The new 
concept of "inherent dynamics" is introduced, and the results from applying 
this concept to the model described in chapter 3 are discussed. 

Chapter 5 deals with the symmetric hopping model. A new numer- 
ical method for investigating the model in the extreme disorder limit is 
developed. Numerical results are presented and compared with 3 analytical 
approximations . 

The three main chapter (3-5) contains individual conclusions. 
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2.3 Papers 

The following papers contain results obtained for the two models discussed 
in the present thesis: 



• Paper I: Effective one-dimensionality of universal ac hopping conduc- 
tion in the extreme disorder limit. J.C. Dyre and T.B. Schr0der. Phys. 
Rev. B. 54 14884 (1996). 

• Paper II: Hopping in a supercooled binary Lennard-Jones liquid. T.B. 
Schr0der and J.C. Dyre. J. Non-Cryst. Solids. 235-237 331 (1998). 

• Paper III: Crossover to potential energy landscape dominated dynamics 
in a model glass-forming liquid . T.B. Schr0der, S. Sastry, J.C. Dyre 
and S.C. Glotzer. J. Chem. Phys. 112 (2000) in press. 

(http: \\xxx.lanl.gov: cond-mat /9901271) 

• Paper IV: Scaling and Universality of ac Conduction in Disordered 
Solids. T.B. Schr0der and J.C. Dyre. Phys. Rev. Lett. 84 310 
(2000). 

• Paper V: Universality of ac conduction in disordered solids. J.C. Dyre 
and T.B. Schr0der. Rev. Mod. Phys. 72 (2000) in press. 



Paper I reports numerical results for the symmetric hoping model, which 
was done mainly during my master thesis [32,33]. Paper II-V report results 
obtained as part of the present Ph.D. thesis. Paper IV and V was written 
after the conclusion of my Ph.D. (summer 1999). 



CHAPTER 2. INTRODUCTION 



Chapter 3 

A Model Glass-former 



A glass-forming binary Lennard-Jones liquid is investigated by molecular 
dynamics under (pseudo) equilibrium conditions. In this chapter we report 
'standard' results for this model liquid, i.e. pair-correlation functions, mean 
square displacements, intermediate scattering functions, etc. The main re- 
sult is, that at low temperatures hopping is present in the system as indi- 
cated by a secondary peak in 4nr 2 G s (r, t), where G s (r,t) is the van Hove 
self correlation function. It has not been possible to identify a temperature 
range, where the asymptotic predictions of the ideal mode coupling theory 
are fulfilled. Some of the results reported in this chapter are contained in 
paper II and paper III. 



3.1 Model and Method 



The system investigated in the present work, has earlier been shown to 
exhibit hopping by Wahnstrdm l|j. By "hopping" is here meant, that par- 
ticles behave like illustrated in figure 3.1; The particle shown stays relatively 
localized for a period of time, and then move some distance, where it again 
becomes localized. The presence of hopping was the main motivation for 
investigating this system, and this kind of dynamics will be discussed in 
more detail in this and the following chapter. 
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CHAPTER 3. A MODEL GLASS-FORMER 




Figure 3.1: Trajectory (in the x-z plane) of a "hopping" particle. First, the 
particle stays relatively localized for a period of time, seemingly oscillating 
around a position in the left side of the figure. After this, the particle move 
more or less directly to a new position in the right of the figure, where it 
again become localized. The temperature is T = 0.59, and the elapsed time 
is At = 160 (see later for units of these numbers). 



In the work of Wahnstrom, the system was investigated both above and 
below the glass-temperature, T g , i.e. the system was allowed to fall out 
off equilibrium as it was cooled. As a consequence, it was unclear whether 
the hopping seen was a feature of the equilibrium liquid, or if it was a 
consequence of non-equilibrium dynamics. Here we attempt to keep the 
system equilibrated at all temperatures, i.e. the equilibration time is made 
longer and longer as the temperature is lowered. 

The "Wahnstrom system" is a binary mixture of N = 500 particles, with 
50% particles of type A, and 50% particles of type B[]. The particles interact 
via the pair-wise Lennard-Jones potential, where the parameters depend on 
the types of the two particles involved (a and (5): 



V aP {r) = Ae aP [[^) -(^j j (3.1) 
The forces are given by the negative gradient of the potential: 

F a/3 (r) = -VV aP = -^ T - (3.2) 



48e a/3 //0- a /3\ 14 1 f<?af3\ 8 



(^a/3) 2 \\ r J 2 



(3.3) 



The exact number of particles used in the simulations are: Na = 251, and Nb = 249 



3.1. MODEL AND METHOD 
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V a p(r) is characterized by a minimum at V a p(2 l l & o a p) = — e a p, steep 
repulsion at shorter distances, and a weaker attraction at longer distances. 

The length-parameters used in the Wahnstrom system areQ: cfaa = 1, 
&bb = 1/1-2 « 0.833 and oab = (o'aa + crb)/2. The energy-parameters 
are all identical: 6aa = ^AB = ^bb = 1- The masses of the particles are 
given by tjia = 2, and tub = 1- The length of the sample was L = 7.28, 
which gives a (reduced) density of p = N/L 3 = 1.296. Times are reported in 
units of r = {mB<y\ A j '48eA^) 1//2 (This expression contains an error in paper 
II). When comparing simulations of a model-liquid like the one described 
here with experimental results for real liquids, it is customary to use "Argon 
units", i.e. parameters used when modeling Argon atoms by the Lennard- 
Jones potential; a = 3.4A, e = 120Kk B , and r = 3 x 10~ 13 sec Jl7| ]. 

The potential is "cut and shifted" |l[] at r = 2.5a a ^, which means that 
it is set to zero for r > 2.5a a /3, and \V a p(2.5o a p)\ « 0.016 is added to the 
potential for r < 2.5a a ^. This makes the potential continuous at r = 2.5cr Q , / g, 
while the force is not. The equations of motion are integrated with periodic 
boundary conditions using the Leap-Frog algorithm ]l]] (which is simply a 
discretizaition of Newton's second law) with a time step of O.Olr. 

Three independent samples were used, each initiated by generating a 
random configuration followed by equilibration at a high temperature (T = 
5.0). The cooling was done by controlling the total energy of the system, 
which was done by scaling the velocities. An equilibration run of length t eq 
was then performed, followed by a production run of the same length. If the 
samples was determined not to be equilibrated properly, t eq was doubled, 
by "degrading" the production run to be part of the equilibration run, and 
making a new production run twice as long. This procedure was continued 
until the samples was determined to be equilibrated. The criteria used for 
determining if the samples were equilibrated will be discussed later. Note, 
that by this procedure, it is the total energy (E to t = E pot + E^ in ) that is 
controlled. The reported temperatures and pressures are computed as time- 
averages over the instantaneous temperature and pressure respectively 

T(t) ^ (3.4) 
P{t) = P T(t) + W(t)/V, ^(tl^^r^F^ (3.5) 

i j>n 

where W(t) is the virial and the summing is over all pair of particles. 



2 We follow |34| and term the large particles "A", and small particles "B" 
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CHAPTER 3. A MODEL GLASS-FORMER 



Some of the results for the system described above, will be compared with 
results from a different binary Lennard- Jones mixture (the "Kob & Andersen 
system"), which has been investigated in a number of papers [17, 18,^4-36]. 
This system is a 80:20 mixture, with gaa = 1, &AB = 0.8, obb = 0.88, 
gab = 0.8, caa = 1-0, €ab = 1-5, e B B = 0.5, and m A = m B = 1. 



3.2 Static Results 

Figure |3.2j shows as a function of temperature a) the potential and total 
energy per particle, and b) the pressure vs. Temperature. The error-bars are 
estimated from deviations between the three independent samples, which are 
found to be in reasonable agreement. This is a necessary (but not sufficient) 
condition for the system(s) to be equilibrated. Note that its the total energy 
which was set to be equal for the three samples, and consequently there are 
small deviations in the measured temperatures. 



-3.0 




5 I i i i i i i i i i i i i i 

' 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 
Temperature 



Figure 3.2: a) Total and potential energy per particle as a function of tem- 
perature, b) Pressure vs. temperature. In both a) and b) (as in the following 
if nothing else is mentioned) error-bars are estimated from deviation between 
three independent samples. 



3.2. STATIC RESULTS 
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The specific heat capacity at constant volume, C v , is given by p[ 



1 dEtot , v 

a = iv^r" ( 3 - 6 ) 



From figure 3.2a we see that C v increases as the temperature is lowered. 
Note that there is no indication of a glass-transition, since this would be 
indicated by a sharp decrease in C v when cooling below the glass transition 
temperature, T g , as seen in the results by Wahnstrohm. The heat capacity 
at constant volume will be discussed further in section ^. 



In figure 3^ is shown the pair-correlation functions, 9aa{t), 9ab(?~), and 
9bb i r ) > at three of the temperatures simulated. The pair-correlation func- 
tion, g a p{r), is the average relative density of particles of type (3 around 
particles of type a []37j. It is normalized to be 1 at large distances, r, were 
the correlations disappears. At the highest temperature (T = 1.06) all 
three pair-correlation functions is seen to look like 'typical' high temper- 
ature pair-correlation functions, with a sharp first neighbor peak, a more 
rounded second neighbor peak, etc. As the temperature is lowered the first 
neighbor peak becomes sharper and the second peak splits into two peaks. 
The splitting of the second peak upon cooling is often seen in super-cooled 
liquids (see eg. [p|, |17[ ) . 

The parameters used in the potential does not energeticly favor the mix- 
ture of A and B particles (e^B is not larger than eaa and £bb), which might 
cause the system to phase separate at sufficiently low temperatures. In the 
event of phase separation, we would expect the area of the first peak in 
9ab{t) to decrease at the temperature where phase separation starts occur- 
ring. There is no indication of this in figure 



In figure 3A is shown the static structure factors, Saa(q), Sab(q), and 
Sbb(q), corresponding to the data shown in figure |3.3| . The static structure 
factor is the (3 dimensional) Fourier transform of pair correlation function, 
and thus contains the same information. However the finite sample size 
introduces some features in S a p(q) which are clearly unphysical (see eg. |jj| ). 
By recalculating S a p{q) from g a p{r) with a cut-off in r that is less than L/2 
one can get an idea about which features of S a p{q) are due to finite-size 
effects. Doing this shows, that at q- values to the left of the first peak (q as 6, 
depending on type of correlation) S a /3(q) is dominated by finite-size effects, 
which is why most of it is not shown. Also the small "wiggles" seen at q ~ 10 
is a finite-size effect. The splitting of the peaks at q ~ 15 is however not a 
finite-size effect. 
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Figure 3.3: Pair correlation function at a high (T = 1.06), medium (T = 
0.66), and low temperature (T = 0.59); a) A-A correlation, gAA( r ), b) A-B 
correlation, gAB(r), c) B-B correlation, gBB{f) 




Figure 3.4: Static structure factor at a high (T = 1.06), medium (T = 
0.66), and low temperature (T = 0.59); a) A-A correlation, Saa(q), b) A-B 
correlation, Sab(q), c) B-B correlation, Sbb(q) 



18 CHAPTER 3. A MODEL GLASS-FORMER 

3.3 Mean Square Displacement 



In figure |3,5| we show the mean square displacement of the a) A particles, 
{r 2 (t))A, and b) B particles, (r 2 (t))#. At short times (t < 1) a ballistic 
regime ((r 2 (t)) Q oc t 2 ) is seen for both A and B particles, at all temperatures. 
The ballistic regime is simply a consequence of the velocities being constant 
at these short time scales: 



(r 2 (t)) a = ((Vt)\ 



(V 2 



(3.7) 



In the insets in figure ^5| is shown (V 2 ) a calculated by eq. 
also be calculated directly from the temperature: 



]. (V 2 ) a can 



2E kin Nm a ((V/T) 2 ) a 48m Q (y 2 ) ( 



3(iV-l) 



3(N- 1) 



3(1 - l/N) 



(3.8) 



In the insets in figure 3J3 is also shown (V 2 ) a as calculated from eq. |3.8| . 
The excellent agreement between eq. |3.7| and 3.8 (with no fitting involved) 



is not a big surprise, but it gives confidence in the argument leading to eq. 
[Q, and acts as a consistency check. 



In figure 3.5 a diffusive regime ((r 2 (t)) a oc t) is seen at long times, for 
all temperatures. As the temperature is lowered the time scale at which 
the diffusive regime sets in increases, and a plateau is seen to evolve be- 
tween the ballistic and diffusive regimes. This behavior is typical for what 
is seen in simulations of super-cooled liquids, and the plateau is argued to be 
associated with particles being trapped in local "cages" consisting of their 



neighbors [17 



The vertical dashed lines in figure 3.5 identifies the time t±, defined by 
(r 2 {ti)) a = 1. The significance of this time will be discussed in connection 
with figure |3.6|. 





Figure 3.5: The mean square displacement of a) the large particles, (r 2 (t))A, 
and b) the small particles, (r 2 (t))#. Insets show (V 2 ) a calculated from eq. 
|3.7| and |3.8| , and demonstrates the agreement between the two equations. 
The vertical dashed lines identifies {r 2 {t\)) a = 1. 



20 



CHAPTER 3. A MODEL GLASS-FORMER 



Fitting (r 2 (t)) a in the diffusive regime to the Einstein relation: 



(r 2 (t)) a = 6Dt 



(3.9) 



we determine the diffusion coefficient, D, as a function of temperature (and 
particle type). We postpone the discussion of the temperature dependence 
of D to section |3.6|, where it will be discussed together with other measures 



of the time scales involved (eg. ti). Following Kob and Andersen [17| we in 
figure 3.6 present the data for (r 2 (t)) a as a function of 6Dt. In the diffusive 



regime the data for all temperatures should, according to eq. |3.9| , fall on 
the line 6Dt, which is seen to be the case for (r 2 (t)) a > 1. In other words 
the time t\, defined by {r 2 (ti)) a = 1, can be viewed as marking the onset 
of the diffusive regime. The fact that the diffusive regime is reached at all 
temperatures (compare fig. |3.5| ), is a necessary (but not sufficient) condition, 
for the system(s) to be in equilibrium. 

As is the case for the data presented by Kob and Andersen, the data 
in figure seems to indicate a universal behavior; As the temperature is 
decreased (r 2 (t)) a follows the thick dashed curve in figure 3J3 to lower and 



lower scaled times. The thick dashed curves are fits to the fitting function 
used by Kob and Andersen |l7j (see table EO|): 



+ A(Dt) b + 6Dt 



(3.10) 



The parameter b is the so-called "von-Sweidler" exponent, which is related 
to the dynamics at intermediate times, i.e. between the plateau and the 
diffusive regime. According to the asymptotic predictions of the ideal mode 
coupling theory (MCT) [||,||, the von-Sweidler exponent, b, should be 
independent of the particle type. Kob and Andersen argue that the two 
values they find, 0.48 for the A particles and 0.43 for the B particles, are 
within reasonable agreement. This is not the case for the b- values found 
here; 0.17 ± 0.03 for the A particles and 0.27 ± 0.01 for the B particles. 



Type 


r c 


A 


b 


A 


0.131 ±0.013 


0.071 ±0.002 


0.17 ±0.03 


B 


0.166 ± 0.003 


0.096 ± 0.003 


0.27 ±0.01 



Table 3.1: Parameters found by fitting equation |3.10j to the "universal" 

o confidence 



curves in figure 3J3 (thick dashed lines). Error-bars are 
intervals, as reported by the fitting routine in Gnuplot. 



3.3. MEAN SQUARE DISPLACEMENT 21 




Figure 3.6: The mean square displacement plotted as function of 6Dt (same 
data as in figure For (r 2 (t)) a > 1 the data for all temperatures fall 

on the straight dotted line marking diffusive behavior (eq. |3.9| ), The thick 
dashed lines are fits to eq. [HO] 
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3.4 The van Hove Correlation Function 

To characterize the dynamics in more detail, we calculate the van Hove self 



correlation function [37| (where the sum is over particles of type a): 



which is the probability density of a particle of type a S {^4, B} being 
displaced by the vector r, during the time interval t. In an isotropic sys- 
tem G sa (r, t) does not depend on the direction of r, and the probability 
distribution of a particles being displaced the distance r is A-Kr 2 G sa {r, t). 
The mean square displacement, (r 2 (t)) a , is given by the second moment of 
4:TTr 2 G sa (r, t), and the later thus gives a more detailed view of the dynamics. 

In figure ^S7j\ is plotted 47tr 2 G s A{r, t\) and 4-7rr 2 G s s(r, t\), i.e. the dis- 
tribution of distances moved by A and B particles respectively in the time 
interval t\ , which in the previous section was demonstrated to be at the onset 
of the diffusive regime. The thick curve is the Gaussian approximation: 

/ 3 \ 3/2 ( 3r 2 
G s{r,t) = - . 2 , .. exp 



2vr(r 2 (t)) J 1 \ 2{r 2 {t)) 

which is seen to be reasonably fulfilled at high temperatures. As the temper- 
ature is lowered the results starts graduately deviating from the Gaussian 
approximation, and a shoulder builds at the average inter particle distance, 
r ~ 1.0, which at T=0.59 becomes a well-defined second peak. The second 
peak, observed also in other model liquids, is interpreted j|,|,13-[l7j as sin- 



gle particle hopping (see figure |3.1[) , i.e. the particles stay localized in their 
local cages a certain time (first peak), then moves more or less directly out 
to distances approximately equal to the inter particle distance (second peak) 
(This behavior was also reported in |]l2)). To what extent the hopping be- 
havior is related to activated hopping over energy barriers will be discussed 
in chapter [|. 

Note that when analyzing the dynamics in more detail by means of the 
van Hove self correlation function (fig. ^^), the universality seen in the 
mean square displacement (fig. |3.6|) is completely gone; The way the system 
"achieves" the diffusive behavior (i.e. the dynamics at t\) changes qualita- 
tively from high temperature (Gaussian) to low temperature (Hopping). 
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Figure 3.7: The distribution of particle displacements, A7rr 2 G sa (r, ti), 
where t\ is defined by (r 2 (t\)} a = 1, see figure |3.5| . At high temperatures 
the Gaussian approximation (thick curve) is reasonable fulfilled, whereas at 
the lowest temperature a secondary peak is present, indicating that hopping 
is present in the system (see text). The small particles (type B) are seen to 
have a larger tendency to exhibit hopping, compared to the large particles. 
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Figure 3.8: The non-Gaussian parameters, a.2A{t) and 02b(*)- The dashed 
line is a fit to a power-law in the "universal" regime before the maximum. 

The deviation from Gaussian behavior is often analyzed in terms of the 
Non-Gaussian parameter ]T7|,|3T 



a 2 (t) 



<r 4 (i)) - (rHt))Gauss 



<r 4 (i)) 



Gauss 



3(r 4 (t)) 
5(r 2 (t)) 2 



1 



(3.12) 



which is the relative deviation between the measured (r 4 (t)} and what it 
would be for a given value of (r 2 (t)), if G s (r, t) was Gaussian: (r 4 (t)) Gauss = 
5{r 2 (t)} 2 /3. In figure [3.8| is shown a 2 for the A and B particles respectively. 
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Like in [17,38] we see a universal behavior in the way a*2A {t) increases, 
before it reaches its maximum. Following |}8[ a power-law fit is done in this 
regime, from which is found exponents of 0.64 and 0.67 for the A and B 
particles respectively. The exponent found in |38| is 0.4. As noted in [0] 
the universality seen in a2(t) is different from the one found in the behavior 
of (r 2 (t)) a (fig. |3.6| ); The universality seen in a^{t) does not involve any 
scaling. 

The time t* which is defined as the position of the maximum in 0:2 (i), can 
be interpreted as the time where the dynamics deviates most from Gaussian 
behavior [^,|4^]. In figure |3.9| is shown for three temperatures, the time 
development of 4nr 2 G s B{r,t). At each temperature 47rr 2 G s s(r, t) is shown 
at a time close to[] t* (indicated by the arrows), and 2,4,8, and 16 times that 
time. At T = 1.06 the distribution of particle displacements is seen to be 
characterized by a single peak, which "spreads out" without any qualitative 
change in the shape, as time increases. This is the typical behavior for 
liquids at high temperatures |37j. At T = 0.66 the behavior is seen to be 
qualitatively similar, except for a weak indication of a shoulder at r ~ 1.0. 
The time development of 47rr 2 G s s(r, t) at T = 0.59 is qualitatively different 
from what is found at higher temperatures; The first peak decreases while 
the position is almost constant, and at the same time as the second peak 
starts building up. This supports the "hopping" interpretation of the second 
peak given above; The particles escape their local "cages" (the first peak), by 
"suddenly appearing" at approximately the inter particle distance (second 
peak). After the last time shown here, the second peak starts decreasing. 
This is interpreted as a consequence of particles hopping again. 

With the regards to the interpretation of t* mentioned above, one should 
note that for T = 0.59 the second peak in 47rr 2 G< j s(r, t) builds up after the 
time t* , i.e. after 0:2 (i) has its maximum. A more detailed analysis of the 
non-Gaussian behavior should include the higher order analogues of 0:2 (t) 
(involving higher moments of 47rr 2 G SQ (r, t)) p9| , but this approach will not 
be pursued further here. 



3 The reason that t* itself is not used here, has to do with at what time-steps configu- 
rations are stored, and thus which time differences are easy accessible 
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Figure 3.9: Time development of the distribution of particle displacements 
for the B particles, 47rr 2 G s B(r,t), at the same temperatures as in figure |3~lj| . 
At each temperature 4irr 2 G s B(r,t) is shown at a time close to t* (smallest 
time indicated by arrows), and 2,4,8, and 16 times that time. 
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The van Hove distinctive correlation function, Goapir, t), is the time 
dependent generalization of the pair correlation function Jl3j : 

G DaP {r,t) = -——^^{Sfoit) - n(0) - r)) , (3.13) 
y/N a N p , =1 . =1 

where the first sum is over particles of type a, the second sum is over particles 
of type (3, and % ^ j if a = (3 (this last term is contained in G sa (r,t), see 



equation 3.11 ). GD a p{r,t) is the relative density at time t of /3-particles at 
r, given that there at time t = was a a-particle at r = (normalised to 
be 1 if there is no correlation). 



In figure 3.10| is shown GoBB( r ,t) for the same temperatures and times as 



in figure Included in figure |3~1C| is also Gdbb{t, t = 0) = gBB{f)- Except 
for the lowest values of r, the features of the pair correlation function is seen 
to approach the long-time value Gu a p{r, t) = 1 in a smooth manner, which is 
the typical high-temperature behavior |37[| . As the temperature is decreased, 
a significant "overshoot" develops at the small r-values^j, indicating that 
there is an excess probability that when a particle jumps, it does so to a 
position previously occupied by another particle. Thus, while G s B(r,t) (fig. 
|3.9| a) shows that particles has a tendency to jump a distance approximately 
equal to the inter-particle distance, GDap{r,t) give the further information, 
that they have a tendency to jump to a position previously occupied by 
another particle. Similar (but less pronounced behavior) was found in |13|| . 

In a scenario where the dynamics is dominated by particles jumping 
to positions previously occupied by other particles, one should expect that 
moving particles make up correlated string-like objects in the liquid. This 
kind of behavior is found by Muranaki and Hiwatari in a soft sphere system 



[44 1, and by Donati et. al. in the Kob & Andersen system This type 



of behavior will be discussed further in chapter 



4 the data for r < 0.1 is removed, since the statistical noise dominates the data in this 
range 
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Figure 3.10: Time development of Gdrb(?", i) for the same temperatures 

e same as 
0) equals 



as in figure 3.9. For each temperature, the times used here are the same as 
(thick dotted curve), where GDBB(r,t 



in figure 3.6, and t 
the pair correlation function, gBB( r )- 
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3.5 The Intermediate Scattering Function 



The self part of the intermediate scattering function, F s (q,t), is the (3 di- 
mensional) Fourier transform of the van Hove self correlation function ||3^ 



*i(q,t) 



G s {r,t)e-^ r dr = (cos q(r 3 (t) - rj (0))) (3.14) 



F s (q,t) is a relaxation function, normalized to be 1 at t = 0, and it goes 
to zero for t — > oo. Figure p. 11 shows the self part of the intermediate scat- 
tering function for the A particles and B particles; F s A(q m axA = 7.5, t) and 
F S B(qmaxB = 8.1, t). q ma xA = 7.5 and q ma xA = 8.1 are here the positions of 
the first peak in the static structure factor for the A- A correlation (Saa(q)), 
and the the B-B correlation (Sbb(q)), respectively (see figure |3.4|) . As the 
temperature is lowered a two-step relaxation develops, which is what is typ- 
ically seerj0 in glass-forming liquids, both in experiments [46| and in simula- 
tions ||,||,|l8|,||,||,||,|3|,||jt^-||]. The initial relaxation is a consequence 
of the particles oscillating in their cages, while the the long-time (alpha) re- 
laxation is a consequence of particles escaping their cages |3g, [|2[^| . The 
alpha relaxation is often found to be well approximated by stretched expo- 
nentials, f(t) = f c exp(-(t/r Q )' 3 ) 



which is also the case 



in figure [37TT], w here fits to stretched exponentials are shown as dashed lines. 
In Fig. |3,12 we show as a function of temperature the three fitting param- 



eters used in figure 3.11; a) relaxation times, r a , b) stretching parameters, 
(3, and c) non-ergodicity parameters f c . As expected the relaxation times, 
T a A and r a B is seen to increase dramatically as the temperature is low- 
ered. The asymptotic prediction of the ideal mode-coupling theory (MCT) 
is t q = To(T — T c )~ 7 [^1, 22]. The divergence of r a at T c predicted by the 
ideal MCT is never seen in practice. This is argued to be consequence of 
relaxation by hopping taking over close to T c [15, 17, 21, 22, 5(J. This means 
that close to T c the power-law is expected to break down and it should not be 
fitted in this regime. Unfortunately we have no independent estimate of T c , 
and furthermore it is not known how close to T c the power-law is expected 
to hold. Since we have already demonstrated, that hopping is present in the 
system at the lowest temperatures, we can expect that these are close to T c , 
which means that we might run into the problem described above. We note 
also, that the power-law is an asymptotic prediction, i.e. it is expected to 



5 In most experiments, and some simulations, this is seen as two peaks (one for 
each relaxation step) in the imaginary part of the generalized susceptibility: x'i^ - 1 ) = 



ujirF(q,u)), where F(q,u>) is the Fourier transform of F(q,t), see eg. 



Jtibility: 



30 



CHAPTER 3. A MODEL GLASS-FORMER 



break down at high temperatures, but again it is not known how far from 
T c this is supposed to happen. 



1.0 
0.8 
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Figure 3.11: a) The self part of the intermediate scattering function for the 
A particles, F s A(q ma xA = 7.5, t). The dashed lines are fits to stretched expo- 
nentials, f(t) = f c exp(— (i/r Q .) /3 ). The fitting was performed for i > 10 for 
the 2 highest temperatures, and for t > 30 for the rest of the temperatures, 
b) F s B(q m axB = 8.1, t), otherwise as above. 
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Figure 3.12: Fitting parameters used in figure |3.11| . a) Relaxation times 
r a A and T a B- The lines are fits to r a oc (T — T c )~ 7 (see text), b) Stretching 
parameters and /3s. c) Non-ergodicity parameters, f c A and / c s. 
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As discussed above, the temperature range (if any) where the asymptotic 
predictions of ideal MCT is supposed to hold, is not known. Consequently 
the following approach for fitting r a = tq(T — T c )~ 7 has been applied; First 
the fit was done for all 8 temperatures (T cut = 0), then by excluding the 
lowest temperature (T cut = 0.60), and then by excluding the two lowest 
temperatures (T cut = 0.63). 

The results of the fitting-procedure described above, is seen in table 
[3.2j . For both the A and B particles the parameters found for T cu t = 0.60 
and T cu t = 0.63 agree within the error-bars, while they disagree with the 
parameters for T cut = 0. Thus one might argue that the parameters found 
for T cut = 0.60 and T cut = 0.63 describes the "true" power-law, while the 
fit found using T cut = deviates as a consequence of fitting too close to T c . 
Thus our best attempt at a power-law fit is the one achieved for T cut = 0.60 
(i.e. by excluding the lowest temperature), which is the fit shown in figure 
|3.12| for the A and B particles respectively. The values estimated in this 
manner for T c and 7 for the A and B particles are identical within the error- 
bars, as predicted by the ideal MCT. The temperature dependence of r a 
will be discussed further in section [0| 



Type 


T C ut 


TO 


T c 


7 


A 


0.00 


4.2 ±0.6 


0.571 ±0.004 


1.71 ±0.10 


A 


0.60 


5.4 ±0.5 


0.592 ± 0.004 


1.41 ±0.07 


A 


0.63 


5.3 ±0.4 


0.586 ± 0.008 


1.46 ±0.09 


B 


0.50 


2.9 ±0.5 


0.571 ± 0.004 


1.77 ±0.11 


B 


0.60 


3.7 ±0.4 


0.592 ±0.005 


1.47 ±0.09 


B 


0.63 


3.6 ±0.4 


0.584 ±0.010 


1.52 ±0.12 



Table 3.2: Parameters found by fitting r a = tq(T — T c ) _7 to all 8 tem- 
peratures (T C ut = 0), excluding the lowest temperature (T cu t = 0.60), and 
by excluding the two lowest temperatures (T cu t = 0.63). The error-bars 
indicate the 68.3% confidence interval, as reported by Gnuplot. 



The numerical data for r a (T) does not by itself give strong evidence 
for a dynamical transition at the estimated T c = 0.592 ± 0.004; nothing 
special seems to happen at that temperature. However, the agreement with 
the lowest temperature (T = 0.591 ± 0.002), where we clearly see hopping 
(figure 3.7), gives us confidence that there is a dynamical transition close to 
the estimated T c . 

The stretching parameters, (3a and (5b, are seen in figure |3~T^ b to de- 
crease from values close to 1 (i.e. exponential relaxation) at high temper- 
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atures, to values close 0.5 at the lowest temperatures. The values found 
at high-temperatures are somewhat dependent on the time-interval used for 
the fitting, and the real error-bars are thus bigger than the ones shown 
(which are estimated from deviations between the three independent sam- 
ples, with the fitting done in the same time-intervals). However focusing on 
the medium to low temperatures (T < 0.73), there can be no doubt that 
(3 is increasing as function of temperature. This is in contradiction with 
the asymptotic predictions of the ideal MCT, which predicts (3 to be con- 
stant. The ideal MCT does not directly predict the long-time relaxation 
to be stretched exponentials, but it predicts it to exhibit time-temperature 
super-position (TTSP) [21, 2^]. This mean that the shape of the long-time 
relaxation should be independent of temperature, which for stretched expo- 
nentials means that (3 should be constant. 




Figure 3.13: (a) The self part of the intermediate scattering function for 
the A particles, F s a(q = 7.5, t), scaled to test if the long-time relaxation 
exhibits time-temperature super-position. This is not the case. 



34 



CHAPTER 3. A MODEL GLASS-FORMER 



The question of whether (3 is constant or increases towards 1 has led 
to controversy about what is the right fitting procedure |50|,|5l|. Another 
way of checking for TTSP, is by scaling the data appropriately and look 
for approach to an universal curve. This is done in figure 3.13; where 
FsA(lmax,t)/ fcA is plotted versus t/r a A which should be identical for any 
temperature range where TTSP holds. No such temperature range can be 
identified. One might argue, that the scaling approach used in figure 3.13j 
relies on the values of f c estimated from the fits to stretched exponentials, 
which at the high temperatures is dependent on the time-range chosen for 
the fitting. Scaling F sA (q ma x,t) to agree at F sA (q m ax,t) = e" 1 (i.e. without 
dividing by / c ), like done in |l8| does not change the conclusion above; there 
is no indication of TTSP. 

The fact that F sa (q,t) decays to zero (figure [3.11 ), and that the pa- 
rameters describing the alpha relaxation is in reasonably agreement for the 
three independent samples (as illustrated by the error-bars in figure 3.12 ), 
are necessary conditions for the liquid to be in equilibrium. Of all tests for 
equilibrium applied, this was the one that was found to be most sensitive, 
i.e. requiring the longest equilibration times. 

In figure 3.14| is plotted F s A(q,t) at T = 0.59 for q- values 1, 2, 20, and 
q = q-max = 7.5 (thick dashed line). Also shown in figure 3.14 are fits to 
stretched exponentials (dashed lines). The fitting parameters used in 3.14 



is plotted in figure |3.15| . At the lowest q- values the fits are not perfect, but 
except for that the fits are reasonable. Note that in the time-range were the 
second peak builds up in 4nr 2 G sA (r, t), i.e. from t* « 6 • 10 2 to h rj 7 • 10 3 , 
figure p. 14 shows no clear indication of this. S ince F sct (q,t^) is the Fourier 
transform of G sa (r, t), and thus contains the same information, it is possible 
to extract the information about the hopping from F sa (q,t). However since 
we already have an excellent indication of the hopping in 4Trr 2 G sa (r, t), and 
Go a p(r,t), this will not be pursued further here. 
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Figure 3.14: The self part of the intermediate scattering function for the A 
particles, F s A(q,t), at T = 0.59 for g-values 1,2,..., 20, and q = q ma x = 7.5 
(thick dashed line). Thin dashed lines are fits to stretched exponential, 
where the fitting was done for t > 30. 




Figure 3.15: q-dependence of parameters found by fitting stretched expo- 
nentials to F s A(q,t), at T = 0.59, see figure 3.14 
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3.6 Time Scales 



In this section, we compare the measures of time-scales found in the previous 
sections. In figure |37l6| is plotted (6D) -1 (squares, from figure pT5| ) , t\ (plus, 
from figure 3.5), t a (diamonds, from figure |3~12] ) , and t* (circles, from figure 
|3.8] ). If t\ is in the diffusive regime, the relation (6D) _1 = t\ should hold, see 
equation 3JJ . This is seen in figure 3,16 to hold to a good approximation; the 
largest deviation is 8%, with t\ being smaller than (6-D) -1 , thus confirming 
that t\ is a good estimate of the onset of the diffusive regime, as argued in 
section |3~3| . t* is found to be roughly a factor 10 smaller than ti, and (6.D) -1 
(the ratio changes from roughly 6 at high temperatures, to roughly 11 at 
low temperatures). Also shown in figure 3.16| are attempts to fit (6-D) -1 to a 
power-law. According to the asymptotic predictions of ideal MCT, (6-D) -1 
should have the same temperature dependence as T a (apart from a constant 
factor), which means that we should find the same T c and 7. Consequently 
the fitting was done in the temperature range which gave the "best" power- 
law fit to T a , i.e. by excluding the lowest temperature. This fit is shown, 
for the A and B particles respectively, as the full curves in figure 3.16 , and 
the parameters and error-bars are given in table [T^. The fit to the data 
are reasonable, but the estimated T c deviates from the one found from r a , 
which is in contradiction with MCT. The dashed curves in figure 3. IE are 
the results of power-law fits where T c was set to have the value found for 
r a , T c = 0.592. This also fits the data reasonable, but now the exponents 
7 are different from the ones found from r a , thus illustrating that r Q and 
(6-D)" 1 has different temperature dependence (as can be seen directly in 
figure |3~1~6|) . Also Kob and Ander sen find different temperature dependence 
for the diffusion and the relaxation 17, 33]. 



Type 


Tcut 


ro 


T c 


7 


A 


0.60 


43 ±3 


0.574 ± 0.005 


1.40 ±0.09 


A 


0.60 


53 ±4 


0.592 


1.18 ±0.03 


B 


0.60 


34 ± 2 


0.573 ±0.003 


1.33 ±0.05 


B 


0.60 


42 ±3 


0.592 


1.11 ±0.03 



Table 3.3: Parameters found by fitting {QD)~ l = tq(T — T c )~ 7 , excluding 
the lowest temperature (T cut = 0.60). The second set of parameters for 
each type of particles is for T c being forced to have the value found for 
T a , T c = 0.592. The error-bars indicate the 68.3% confidence interval, as 
reported by Gnuplot. 
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Figure 3.16: Time-scales vs. temperature; (6D) _1 (squares, from figure 
|3.5|) , t\ (plus, from figure |3^), t a (diamonds, from figure |3.12| ), and t* 
(circles, from figure |3.8| ). The lines are power-law fits to (6-D) -1 (see text). 
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Figure 3.17: Arrhenius plot of the relaxation times, T a . Full curves: fits 



to the Vogel-Fulcher law (see table 3^4), extrapolated to the laboratorie 



glass-transition (see text). Dashed curves: power-law fits from figure 3.12. 



Figure |3. 17 shows the alpha relaxation times T a A and r a B , in the so-called 
"Arrhenius-plot" , i.e. as a function of 1/T and with a logarithmic y-axis. In 
this plot, a Arrhenius dependence of the relaxation time, r a = Toex.p(E/T) 
would give a straight line, which is clearly not the case. Or in the terminology 
of Angel ]52fl , the liquid is "fragile" (non- Arrhenius) , as opposed to "strong" 
(Arrhenius) . Included as full lines in figure [3.17| are fits to the Vogel-Fulcher 
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law, r a = ro exp(^4/(T — To)), which is often seen to fit the relaxation times 
over a range of temperatures in experiments. The fits are extrapolated, 
to include the time r a ~ 3 • 10 15 which in Argon units (see section |3.1| ) 
corresponds to r a ~ 10 3 seconds, i.e. including the time scales involved 
in the laboratory glass-transition (see introduction). This extrapolation 
over 12 decades in time (from information covering 2.5 decades), should 
of course not be taken too seriously, but is included here to illustrate the 



large differences in time scales. Included as dashed lines in figure 3.17 are 



the power-law fits from figure 3.12 



Type 


TO 


A 


T 


A 


5.2 ±0.9 


0.68 ± 0.07 


0.486 ± 0.009 


B 


3.5 ±0.7 


0.70 ± 0.08 


0.486 ±0.010 



Table 3.4: Parameters found by fitting the relaxation times, r a A and r a B to 
the Vogel-Fulcher law, r a = roexp(^4/(T — To)). 
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3.7 Finite Size Effects 



To estimate finite size-effects we in this section present results from a sample 
with N = 1000 particles which is otherwise equivalent with the samples 
described in the previous sections (i.e. same density and total energy per 
particle). 

Figure |3.18| shows the self part of the intermediate scattering function 
for the A particles for the N=1000 system, together with fits to stretched 
exponential (dashed lines), using the same fitting procedure as for the N=500 
samples (figure |3.11| ). The resulting fitting parameters are shown in figure 
|3.19| together with those found for N=500 (figure 3,12| ). The results for the 
two system sizes are seen to be in reasonable agreement, indicating that 
the results in the previous sections do not depend strongly on system size. 
In IH a decrease in r a by a factor of « 30 was found at low temperatures 
for a soft sphere system by going from N = 108 to N = 10000. At the 
moment it is unclear if this much larger change in relaxation time is due to 
the larger difference in system size, starting from a smaller system, or other 
differences between the two systems. 




10° 10 1 10 2 10 3 10 4 10 5 



Figure 3.18: a) The self part of the intermediate scattering function for 
the A particles, F s A(q m axA = 7.5, i) for a sample with N=1000 particles 



(compare figure 3.11| , T = 0.62 (second lowest) and T = 0.66 (fourth lowest) 
are missing). 
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Figure 3.19: Fitting parameters used in figure 3.1S| for the N=1000 sam- 
ple, a) Relaxation times t q a- Power-law reproduced from fig. 3.11a. b) 
Stretching parameters (3a- c) Non-ergodicity parameters, f c A- 



42 



CHAPTER 3. A MODEL GLASS-FORMER 



3.8 Conclusions 

The original motivation for investigating the binary Lennard-Jones system 
described here, was that it was already demonstrated to exhibit hopping [15]. 



The first goal of the present work was to test if the hopping was still there if 
the cooling was done under (pseudo-) equilibrium conditions. Unfortunately 
there is not a single condition that is known to be sufficient for the system to 
be equilibrated. The different possibilities include: No drift in static prop- 
erties (temperature, potential energy, pressure, etc.), long time dynamics 
being diffusive, relaxation functions such as F s (q,t) decaying to zero. These 
are all necessary conditions, but none of them are (known to be) sufficient. 
In the present work, it was found that the most sensitive condition (i.e. re- 
quiring the longest equilibration time), was that F s (q,t) decays to zero, and 
that it does so in the same manner for the three independent samples, i.e. 
with the relaxation time and stretching exponent being within reasonable 
agreement. Assuming that this condition is sufficient, it is concluded that 
the liquid is in equilibrium at all the temperatures presented here, and thus 
that the hopping found is a feature of the equilibrium liquid. 

Although it was not a goal in itself to test the ideal mode coupling the- 
ory (MCT), the results of the simulations was compared to the asymptotic 
predictions of ideal MCT. This was done for two reasons; I) It provides a 
convenient way of comparing results with other simulations (e.g. Kob and 



Andersen II) It is "common" belief that ideal MCT breaks 

down when hopping dynamics takes over, and the critical mode coupling 
temperature T c thus constitutes an estimate of when hopping should start 
to dominate the dynamics. At first hand the last point seems to work fine; 
the best attempt at a power-law fit to T a gives T c = 0.592 ± 0.004, which 
is close to the lowest of the temperatures simulated, where we clearly see 



hopping (figure |3.7[ ). The fact that hopping also is present at the second 
lowest temperature, might then be taken as an indication that dynamical 
transition from the "MCT regime" to the "hopping regime" is not a sharp 
transition, but takes place over some temperature interval. However, the big 
problem with this scenario is the failure to identify any temperature range, 
where the asymptotic predictions of ideal MCT holds; There is no indication 
of TTSP, and the temperature dependence of the time-scales for diffusion 
and relaxation are clearly different. 

In the full version of ideal MCT (as opposed to the asymptotic predic- 
tions used in the present work), the exponents of the power-laws discussed 
here4 are not free parameters, but are determined from the static structure 



factor [21, 22,54]. Calculating the exponents from the static structure fac- 
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tor is however a difficult numerical task, which has only been done for a 
few systems [54|. This has not been attempted for the system used here. 
In an attempt to include the hopping dynamics in MCT, the so-called ex- 
tended mode coupling theory has been developed introducing an "hopping- 
parameter", as a extra fitting parameter p^] . It has not been attempted to 
apply extended MCT to the system investigated here. 



CHAPTER 3. A MODEL GLASS-FORMER 



Chapter 4 

Inherent Dynamics 



The dynamics of the model glass-forming liquid described in the previous 
chapter, is here analyzed in terms of its "inherent structures", i.e. local 
minima in the potential energy. In particular, we compare the self part of 
the intermediate scattering function, F s (q,t), with its inherent counterpart 
Fg(q, t) calculated on a time series of inherent structures. Fg(q, t) is defined 
as F s (q,t) except that the particle coordinates at time t, are substituted 
with the particle coordinates in the corresponding inherent structure, found 
by quenching the equilibrium configuration at time t. We find that the long 
time relaxation of Fj (q, t) can be fitted to stretched exponentials, as is the 
case for F s (q, t). Comparing the fitting parameters from F s (q, t) and F/ (q, t) 
we conclude, that below a transition temperature, T x , the dynamics of the 
system can be separated into thermal vibrations around inherent structures 
and transitions between these. 

The main conclusions of this chapter can be found in paper III. In the 
present chapter we introduce the concepts of "energy landscape" and "inher- 
ent dynamics" , followed by a summary of paper III. Following this the data 
shown in paper II and some additional data is discussed, and a conclusion 
is given. 
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4.1 The Potential Energy Landscape 



The potential energy landscape flfy|2Q ,|55|-|64| of an atomic system (i.e. par- 
ticles without internal degrees of freedom) is simply the potential energy, as 
a function of the 3N particle coordinates. Letting R denote the 3N dimen- 
sional vector describing the state point of the system (i.e. its position in 
the 3N dimensional configuration space), we write the energy landscape as 
U (R), see figure [O]. The behavior of the system can be viewed in terms of 
the state point, R(i), moving on the energy landscape surface, J7(R). This 
surface contains a large number of minima, termed "inherent structures" by 
Stillinger and Weber [^] . The inherent structures are characterized by zero 
gradients in the potential, and they are thus mechanically stable configu- 
rations. The inherent structures are separated by saddle points acting as 
energy barriers. Each inherent structure has a basin of attraction, in which 
a local minimization of the potential energy (a "quench" ) will map the state 
point R(i) to the corresponding inherent structure R 7 (i). 



U(R)t 



MD configuration, R(t) 




Quench 

[Conjugate Gradient] 



Inherent 
Structure, R ! (t) 

[Stillinger & Weber] 



Configuration space, R 
(3N dimensional) 

Figure 4.1: Basic concepts in the potential energy landscape. The energy 



landscape, U(R), contains a large number (oc exp (kN) [55]) of inherent 
structures, which are local minima in U(R). Any configuration R(t) can 
be mapped by a quench (local minimization of the potential energy), into a 
corresponding inherent structure, R^(t). The fact that the 3iV-dimensional 
configuration space is drawn here as 1-dimensional can be misleading; the 
configuration space is 3iV-dimensional, and should be thought of as such 
(difficult as it might be). One might think of the x-axis in this plot as a par- 
ticular direction in configuration space, connecting two inherent structures. 
It should not be thought of as an order-parameter (Readers who find this 
last sentence strange can safely ignore it). 
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So far all we have done is denning the energy landscape. For a system 
like the binary Lennard-Jones mixture investigated in the previous chapter 
we know ?7(R.) exactly; its simply the sum of the pair-potentials (equation 
|3 . l|) F~|- Of course having a well-defined quantity only helps, if it can tell us 
something about what we are interested in, i.e. in this case the dynamics 
of a glass-forming liquid. This is the main question we deal with in this 
chapter. First we note that for the simulations presented in the previous 
chapter, the dynamics is governed by Newton's second law, which using the 
concepts introduced above can be written as: 

-^R(i) = -M- 1 VC/(R(t)) (4.1) 

where M is a (3Nx3N) diagonal matrix, with the appropriate masses on 
the diagonal. In this sense the dynamics is defined by the energy land- 
scape. This is however obviously not telling us anything new. What we 
want to know is, can we think of the dynamics in terms of the following 
scenario; the dynamics of the liquid is separated into (thermal) vibrations 
around inherent structures and (thermally activated) transitions between 
these. Presented with this scenario, one might very well ask: How can a 
system with constant total energy (e.g. the system simulated in the pre- 
vious chapter) perform thermally activated processes? The answer is that 
one should think of the transitions between inherent structures as involving 
only a few (local) degrees of freedom, while the remaining degrees of freedom 
provides the heath bath. 

In his classic paper from 1969 Goldstein |]l9| argued that there exists 
a transition temperature, which we will term T x , below which the flow of 
viscous liquids is dominated by potential barriers high compared to thermal 
energies, while above T x this is no longer true. Or in other words; below T x 
the dynamics is governed by the "vibrations plus transitions" scenario given 
above. Goldstein gave as a (very) rough estimate, that the shear relaxation 
time at T = T x is on the order of 10 -9 seconds. Later it was noted by 



Angell 1 65], that experimentally it is often found that the shear relaxation 
time is on the order of 10 -9 seconds at the mode coupling temperature, T c . 
This lead to the argument that Goldstein's transition temperature, T x , is 
identical to the mode coupling temperature, T c . A similar argument was 
recently given by Sokolov J66|. 



1 The fact that we know U(R) does not necessarily mean that we can find all the 
inherent structures. This can only be done for small systems (N less than approxemately 
30 depending on other factors such as density and the potential), see eg. [B9|. This is not 
the approach we take here. 
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-6.80 




Figure 4.2: a) Mean value of the quenched energy, (E quenc h) per particle, as 
a function of temperature, b) The heat capacity at constant volume, C V (T), 
and the corresponding "quenched heat capacity", C quenc h(T) (see text). 



Recently, Sastry, DeBenedetti and Stillinger [61] demonstrated that the 
onset of non-exponential relaxation (/? < 1) in a simulated glass-forming 
Lennard-Jones liquid is associated with a change in the system's exploration 
of its potential energy landscape. In figure [4.2| a we present data similar 
to the data presented in J6lj ]. (E quenc h)(T) is here the mean value of the 
energy of inherent structures, found by quenching configurations from a 
normal MD run equilibrated at the temperature T. Like in [61] we find that 



this approaches a constant value at high temperatures, where the relaxation 
becomes exponential {(3 « 1, compare fig. |3.12| ). Similar results was found 
in |67| for the quenched enthalpy (letting the volume change during the 
quench) . 

When plotting the total energy per particle, E tat (T)/N (in figure 3.2), we 



found that the specific heat capacity at constant volume, C V (T), increases 
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as the system is cooled. In figure [4.2| b this is shown explicitly, with C V {T) 
calculated from eq. |3.6| (using central differences). Also shown in figure |4.2f b 
is the "quenched heat capacity", C quenc h(T), calculated in the same way, but 
with (E quenc h)(T) substituted for (E to t)(T). Clearly the increase in C V {T) 
upon cooling is related to the increase in C quenc h{T). 



4.2 The Inherent Dynamics 

The basic idea of the "inherent dynamics" approach is the following (see 
figure |4.3| ); After equilibration at a given temperature, a time series of con- 
figurations, R(t), is produced by a normal MD simulation. Each of the 
configurations in R(t) is now quenched^, to produce the corresponding time 
series of inherent structures, R/(t). We now have two "parallel" time series 
of configurations. The time series R(i) defines the "true dynamics" , which is 
simply the normal (Newtonian) MD dynamics as presented in the previous 
chapter, described by (r 2 (t)}, F s (q,t), etc. In a completely analogous way, 
we define the "inherent dynamics" as the dynamics described by the time 
series JH 1 (t). In other words: If a function describing the true dynamics (eg. 
(r 2 (t)) or F s (q,t)) is calculated by /(R(i)), then the corresponding function 
in the inherent dynamics (eg. (r 2 (i)) 7 or Fg(q,t)) is calculated in exactly 
the same way, except using the time series of inherent structures: /(R / (t)). 

If the true dynamics can be separated into vibrations around inher- 
ent structures, and transitions between these, as stated in the "vibrations 
plus transitions" scenario given above, then the inherent dynamics can be 
thought of as what is left after the thermal vibrations is removed from the 
true dynamics. 



In the bottom part of figure 4.3 the inherent dynamics approach is ap- 



plied to the trajectory of the hopping particle, which was shown in figure 



3.1; All the configurations that were used to plot the true trajectory were 
quenched, and the position of the particle in the resulting time series of 
inherent structures is plotted as the "inherent trajectory". The quenching 
procedure is seen to remove the vibrations in the true trajectory. The mo- 
tion that seems to be left from the vibrations (eg. a jump from x ~ 2.2 to 



x ~ 2.4) will be discussed in section 4.5. 



2 In the present work this minimization was done using the conjugate gradient method 
|68| , which uses a succession of line minimizations in configuration space to minimize the 
potential energy. 
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Figure 4.3: Schematic describing the principle of the "inherent dynam- 
ics" approach. A time series of configurations R(t) of the equilibrated liq- 
uid is generated using normal (Newtonian) MD. Configurations in R(t) are 
quenched to produce their corresponding inherent structures R J (t). Succes- 
sive inherent structures then form a time series which we use to calculate 
"inherent dynamical" quantities such as the inherent mean square displace- 
ment, (r 2 (i)) 7 , and the inherent intermediate scattering function, Fg(q,i). 
Bottom part of the figure: Applying the inherent dynamics approach to the 
trajectory of a hopping particle (see text). 
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Paper II and III explores the possibilities of comparing the true and in- 
herent dynamics of the model glass-former described in the previous chapter. 
The concept of inherent dynamics was first introduced in paper II (without 
using that name). In that paper the true and inherent versions of the mean 
square displacement and the van Hove correlation function, was compared 
on a qualitatively level. Paper III compares the true and inherent version of 
the self part of the intermediate scattering function, which can be done on a 
quantitative level by means of stretched exponentials (see section |3.5| ) . This 
turns out to have important consequences, which is why we here discuss this 
paper first. 



4.3 Paper III 

In this paper the concept of inherent dynamics is applied to the self part of 
the intermediate scattering function, i.e. we compare F s (q max ,t), with its 
inherent counterpart F^ (q max ,t). As explained above F^(q,t) is defined in 
the same way as F s (q,t), except that the normal particle coordinates, ij-(t), 
are substituted by the corresponding coordinates in the inherent structures, 
rj(t) (compare equation |3.14 )P|: 



F/(q,t) = (cosq(rj(t)-ij(O))) (4.2) 



In figure 4b in paper III F^ (q max ,t) is plotted for the A particles, at the 
8 temperatures studied. Here we plot the similar data for the B particles in 



figure 4A In both cases we find that the long time relaxation of i 7 / (q ma x,t) 
can be fitted to stretched exponentials, like is the case for F s (q max ,t). This is 
an important finding, since it enables us to make an quantitative comparison 
between F s (q max ,t) and Fg (q max ,t), by comparing the two sets of fitting 
parameters. In the following, the set of fitting parameters for F s {q max ,t) is 
denoted {r a ,(5,f c }, while the corresponding set for Fj(q max ,t) is denoted 



3 rj(i) is here the 3-dimensional vector describing the position of the j'th particle in 
the inherent structure R, (t). 
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t 

Figure 4.4: The self part of the inherent intermediate scattering function for 
the B particles, F^ B (q = 8.1, t), corresponding to figure 3b in paper III. The 
dotted lines are fits to stretched exponentials: = fl exp(— (V r i) /3 ) 



The key question now is: If the true dynamics follows the "vibrations 
plus transitions" scenario, i.e. if the dynamics can be separated into vi- 
brations around inherent structures, and transitions between these, how is 
{r^,/? 7 ,/^} expected to be related to {r a ,/3, / c }? To answer this question, 
we assume that the initial relaxation in F s (q,t) is due to vibrations, which 
is the widely accepted explanation (see section 3.5). If this is the case, then 
we expect the quenching procedure to remove the initial relaxation (since 
it removes the vibrations), which means that F*(q,i) can be thought of as 
F s (q, t) with the initial relaxation removed. This in turn means, that F*(q,t) 
should be identical to the long time relaxation of F s (q,t), but rescaled to 
start at unity (by definition): {t 7 ,/3 7 ,//} = {r a ,f3, 1}. 

The identity of the relaxation times, and the stretching parameters, 
{t^P 1 } = {r a ,/3}, will probably be true for any "coarse-graining" of the 
dynamics we might apply; if we eg. decide to add a small random displace- 
ment to all the particles (instead of quenching), we would still expect the 
long time relaxation to have the same shape and characteristic time, i.e. 
{t^P 1 } = {t q ,/3}. It is when we find // = 1 we know that the procedure 



4.3. PAPER III 



53 



we have applied removes the vibrations (or in more general terms: removes 
the part of the dynamics responsible for the initial relaxation). 



In figure 5 in paper III the fitting parameters for F s (q max ,t) and (q max ,t) 
are compared for the A particles. Here we show similar data for the B par- 
ticles in figure 4.5. In both cases we find for all temperatures ~ r a . 



The stretching parameters are more difficult to compare at high tempera- 
tures, but they become identical at low temperatures, (3 1 = f3 (within the 
error-bars). Whereas f c is roughly constant as a function of temperature, 
fl is clearly seen to approach unity, as the system is cooled. We have 
thus found evidence, that the the vibrations plus transitions scenario holds 
below a transition temperature T x , as argued by Goldstein. We estimate 
that T x is close to (or just below) the lowest of the temperatures simulated 
(T = 0.591 ± 0.002). 



At the transition we find r Q ~ 3 - 10 3 , which in "Argon units" corresponds 
to T a « 10 -9 seconds (see section 3.1), i.e. Goldstein's estimate of the shear 



relaxation time at T x . This agreement is however probably "to perfect"; 
As mentioned above Goldstein's estimate is very rough, and it regards the 
shear relaxation time and not r a . Consequently an agreement better than 
"orders of magnitude" should probably be considered a coincidence. 



Having found evidence that there exists a transition temperature, T x , 
and estimated the (approximate) position of this, it is natural to proceed 
to check AngelPs proposal, that T x ~ T c . We find that both estimated 
values for T c (0.592 ± 0.005 from relaxation times and 0.574 ± 0.005 from 
diffusion) are in the temperature range where fl is approaching unity. To 
check the proposition that T x m T c on the system used here, is of course 
somewhat problematic, since it doesn't conform very well to the predictions 
of the mode coupling theory, as discussed in the previous chapter. It should 
be noted however, that the arguments given by Angell (and Sokolov), only 
relates to T c as the temperature where power-law fits to experimental data 
tends to break down, i.e. the "usage" of MCT in this argument is similar to 
the way we have estimated T c in the previous chapter, and does not require 
e.g. time-temperature super-position. 



At the end of paper III we present results on the nature of transitions 
between inherent structures. These will be discussed in section 4.5. 
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Figure 4.5: Parameters describing the fits of F s B(q = 8.1, t) and F^ B (q = 
8.1, t) to stretched exponentials; (a) The relaxation time vs. T. The 



solid line is the same power law fit as in figure 3.12 b. (b) The stretching 
parameter (3 1 vs. T. (c) The non-ergodicity parameter, f£ vs. T. 
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4.4 Paper II 

In paper II the concept of inherent dynamics is applied to the mean square 
displacement, (Ar 2 (t)), and the distribution of particle displacements, 
Airr 2 G s (r, t). 

Comparing (Ar 2 (t)) / to (Ar 2 (t)) (figure 1 in paper II), we find that the 
quenching procedure removes the plateau seen in (Ar 2 (t)). This is taken 
as (qualitative) evidence for the "cage" -explanation for the plateau (see 
section |3.3[) . This conclusion is consistent with the conclusion we drew from 
F^(q,t) in the previous section, since the "caging" is what we call vibrations 
in configuration space. 

Comparing 47rr 2 G^(r, t) and 47rr 2 G s (r, t) (figure 2 in paper II), we find 
that the hopping peak is slightly sharper in 4-7rr 2 G^(r, t), and that the first 
peak is moved to the left. A point that is not made in paper II is the 
following; if the particles only moved by either "rattling" in their local cages, 
or hopping approximately an inter-particle distance, then one would expect 
the first peak seen in 4-7rr 2 G s (r, t) to be quenched to a delta-peak at r = 
m 4irr 2 Gl(r,t). Fi gure 2b in paper II shows that this is clearly not the case. 
The reason for this will be discussed in the next section. 



4.5 Transitions between inherent structures 

Having established that our lowest temperature (T=0.59) is close to T x 
(i.e. the dynamics can be thought of as separated into vibrations around 
inherent structures and transitions between these) we are now interested 
in studying the nature of transitions between inherent structures at that 
temperature. We have identified 440 such transitions, by quenching the 
true MD configurations every O.lr (i.e. every 10 MD-steps), and looking for 
signatures of the system making a transition from one inherent structure to 



another. In figure 4.6a is shown the energy of the inherent structures as a 
function of time, E quenc h{t). As expected E 9Uenc f l {t) is found to be constant 
for some time- intervals, and then jump to another level, i.e. the system has 



made a transition from one inherent structure to another. In figure 4.6b is 
plotted as a function of time, the distance in configuration space between 
two successive quenched configurations ]69fl: 



AR T (t) = |R/(t + 0.1) - R/(t)| 



N 



£(r£(t + 0.1)-rj(t) 



(4.3) 
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Figure 4.6: Identifying transition among inherent structures; a) E quenc h vs. 
time, b) AR 1 (equation 4.3 ) vs. time. Transitions between (the basin of 
attraction of) inherent structures is indicated by a jump in E quenc h(t) and 
a corresponding peak in AR ! (t). 



Each jump in E quenc h(t) is associated with a peak in AR J (t), indicating 
the system has moved to a different inherent structure. A transition might 
in principle occur between two inherent structures with exactly the same 
(quenched) energy. Such a transition would not be seen in -£/guenc/i(^)> and 
consequently we use AR 1 ^) to identify transitions. The distribution of 
log 10 (Ai? / (t)) is shown in figure 4.7. The distribution is seen to have two 
separated peaks. The peak to the left (centered around log 10 (Ai? 7 (t)) —3) 
is due to numerical uncertainties; two configurations within the same basin 
of attraction is not quenched to the exactly the same configuration. The 
peak on the right is the one containing the physically interesting transitions, 
which is seen to have AR 1 on the order of unity. In the following we use the 
condition AR 1 > 0.1 to identify the (physically interesting) transitions. 
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Figure 4.7: Distribution of log 10 (AR I (t)). This distribution is not nor- 
malized, i.e. the y-axes indicates how many quenches gave values of 
log 10 (AR I (t)) i n the interval on the x-axis. 



Figure 4.7 shows that most of the quenches results in not finding a 
transition (left peak), and only a small fraction results in actually find- 
ing a transition (right peak). Or in actual numbers: Doing 7500 quenches 
we found 440 transitions. The quenching procedure takes a considerable 
amount of time (corresponding to approximately 1000 MD steps), which is 
why we haven't simply continued this procedure to find a larger number 
of transitions. The data presented in the following is averaged over the 
440 transitions found using the "brute force" method described above, and 
should be considered preliminary. (NOTE: In paper III results from 12000 
transitions are reported. The results are similar to the results presented 
here, except of course with less noise). 

For each transition, we monitor the displacements of the particles from 
one inherent structure to the other. The distribution of all particle displace- 
ments is shown in Fig. 4. £ . While many particles move only a small distance 
(r < 0.2) during the transition, a number of particles move farther, and in 
particular, we find that the distribution for r > 0.2 is to a good approxima- 
tion exponential. At present we have no explanation for this. The dotted 
curve is a fit to a power-law with exponent —5/2, which is a prediction from 
linear elasticity theory jT0| , [7l] , describing the displacements of particles in 
the surroundings of a local "event". This power-law fit does not look very 
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convincing by it self, but we note that the exponent was not treated as a 
fitting parameter (i.e. only the prefactor was fitted), and the power-law is 
expected to break down at small displacements, since these corresponds to 
distances far away from the local event, and is thus not seen in our finite 
sample. From the change in behavior of p(r) at r ~ 0.2, it is reasonably to 
think of particles with displacements larger than 0.2, as those taking part 
in the local event, and the rest of the particles as being in the surroundings, 
adjusting to the local event. Using this definition it is found that on average 
approximately 10 particles participate in an event. 




0.2 0.4 0.6 0.8 1 1.2 1.4 



r 



Figure 4.8: Distribution of particle displacements during transitions be- 
tween consecutive inherent structures at T=0.59. Full curve is a fit to an 
exponential, for r > 0.2. The dotted curve is a fit to a power-law with 
exponent —5/2, as predicted by linear elasticity theory (see text). 



Figure 4.5 has two important consequences with regards to points dis- 
cussed earlier in this thesis; i) The distribution of particle displacements 
during transitions shows no preference for displacement of the average inter- 
particle distance (~ 1 in the used units). This shows that the hopping in- 
dicated by the secondary peak in 4-7rr 2 G s (r, t) (figure 3/7 and |3.9| ) at low 
temperatures is not due to transitions over single energy barriers. A (cor- 
related) sequence of the these is needed, to "build up" the secondary peak. 
This is consistent with the behavior seen in the inherent trajectory in figure 
[4.3j; The jump does not happen in one step, but through a number of "in- 
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termediate" inherent structures, ii) Particles in the surroundings of a local 
event are displaced small distances, to adjust to the larger displacements 
occurring in the local region of the event itself. This kind of motion is very 
hard to detect in the true dynamics, since it is dominated by the thermal 
vibrations. Presumably this kind of motion is the reason why the inherent 
trajectory in figure |4.3| still retains some motion "within" the vibrations; as 
a consequence of an event in the surroundings, the particle starts vibrating 
around a position that is slightly displaced. This view of the dynamics is 
also consistent with the fact, that the first peak in 47rr 2 Gg (r, t) is not a delta 
function in r = (see discussion in section |4.4j ) . 




Figure 4.9: Before (light) and after (dark) one typical transition, all the 
particles which move a distance greater than 0.2. The cooperative, string- 
like nature of the particle motions during the inter-basin transition can be 
clearly seen. 

From visual inspection of a number of the identified transitions it was 
found, that these are cooperative and string-like in nature. By visual in- 
spection is here meant, that the position of particles moving more than 0.2 
during the transition, is plotted before and after the transition, see figure |Q| . 
String-like motion has been found to be an important part of the dynamics 
of glass- forming liquids. It is the natural consequence of particles hopping 
to positions previous occupied by other particles (as concluded from figure 
3.10| ), and it was found (and quantified) in the Kob & Andersen system, 
when looking at the "mobile" particles fl36|] . In paper II strings was also 
found when looking at how particles was displaced (in the inherent struc- 
tures) during the time interval t* (figure 3 in paper II). In paper II this 
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was described as "vacancy hopping"; one particle jumps, leaving room for 
another particle to jump, etc. The finding that also transitions over energy 
barriers are associated with strings, indicates that the vacancy hopping in- 
terpretation might not be correct; it seems to indicate that (at least some 
of) the strings are really cooperative in nature, i.e the particles in the string 
move at the same time. Further (quantitative) investigations are obviously 
needed to answer this question. 

4.6 Conclusion 

In this chapter we have presented results from analyzing the dynamics of a 
model glass-forming liquid in terms of its potential energy landscape. We 
did so by introducing the new concept of "inherent dynamics" , which can be 
thought of as a course-graining of the true dynamics, where the part of the 
dynamics related to vibrations around single inherent structures is removed. 

Comparing the self intermediate scattering function, F s (q,t), with its 
inherent counterpart, Fg(q,t), we found direct numerical evidence for the 
existence of a transition temperature, T x , below which the true dynamics 
is separated into vibrations around inherent structures and transitions be- 
tween these (the "vibrations plus transitions scenario"). We thus confirm 
the "energy landscape" picture, which is (at least) 30 years old. Given the 
fact the energy landscape does exist (since its simply the potential energy as 
function of the particle coordinates) and it does have a number of local min- 
ima (inherent structures), it is not surprising that the dynamics becomes 
dominated by the energy barriers at sufficiently low temperatures. What 
we have done here using the concept of inherent dynamics, is to provide 
direct numerical evidence for this, and we have shown that this regime can 
be reached by (pseudo-) equilibrium molecular dynamics (for the particu- 
lar system investigated here). To our knowledge this is the first time such 
evidence has been presented. 

The fact that we have been able to cool the system, under equilib- 
rium conditions, to temperatures where the separation between vibrations 
around inherent structures and transitions between these is (almost) com- 
plete, means that it now makes sense to study the individual transitions over 
energy barriers, since these in this regime are "significant". There is a lot of 
interesting questions to investigate regarding these transitions, and we have 
here only investigated a few of these. Specifically, we have not determined 
the energy barriers, but only compared the two inherent structures involved 
in a particular transition. This is an obvious point for further investigations. 



Chapter 5 

The Symmetric Hopping 
Model 



The symmetric hopping model is introduced, and three analytical approxi- 
mations for calculating the frequency dependent diffusion coefficient, D(s), 
in the extreme disorder limit (low temperatures) of the model is described; 
the Effective Medium Approximation (EMA) , the Percolation Path Approx- 
imation (PPA), and the Diffusion Cluster Approximation (DC A). DC A is a 
new approximation, developed by my supervisor Jeppe C. Dyre (See paper 
IV and V). Two numerical methods for calculating D(s) in the extreme dis- 
order limit is discussed. The first method is derived from the mean square 
displacement and is equivalent with the method of the ac Miller- Abrahams 
(ACM A) electrical equivalent circuit. The second method (VAC) is derived 
from the velocity auto correlation and is a new method. Numerical results 
using the VAC method are compared to the three analytical approximations, 
and previous results from the ACMA method. 

The main results in this chapter are the development of the VAC method 
(section |5.5[ ), and the numerical results it leads to (section |5]6| and |5.7| ). 
Results in this chapter are published in paper IV and V. 



5.1 The Symmetric Hopping Model 

The symmetric hopping model [23 — 28 1 is defined in the following way: A 



particle 'lives' on the sites of a P-dimensional regular lattice (see figure 5A. 
for a 1-dimensional illustration), where all the sites has the same energy 
(which we set to 0). The particle jumps over energy barriers connecting 
nearest neighbor sites, with jump rates (probability per unit time) given by 
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F(k —>■«) = ro exp(— (3Eki), where Tq is the (constant) "attack- frequency" , 
/3 = {^bT)^ 1 and E^i is the energy barrier between the two sites. The energy 
barriers are chosen randomly, from a probability distribution, p(E), (to be 
specified). Jump rates between sites that are not nearest neighbors are zero. 
It follows from the above, that T(i — > A;) = T(k — > i), i.e. the jump-rates are 
symmetric. In the following we will use T(E) = Tq exp(— (3E), and denote 
the lattice constant a. 




Figure 5.1: The symmetric hopping model in 1 dimension. A particle jumps 
between nearest neighbor sites, by crossing energy barriers connecting these. 
The jump-rates are given by T(i — > i + 1) = Tq exp(— /JE^j+i), where -E^i+i 
is the energy barrier between the two sites. 



If all the energy barriers are identical, i.e. p{E) = 5(E — Eq), the model 
describes diffusion in an ordered structure, and one finds normal diffusive 
behavior: 

(AX 2 (t)) = 2Dt (5.1) 

where (AX 2 (t)) is the mean square displacement along the x-direction, and 
D = cl 2 T(Eq) is the diffusion coefficient pBfl . The average in equation 5.1 



can either be a time-average or a ensemble average. We will in the following 
use the ensemble average, which is characterized by all the sites having the 
same probability (since they have the same energy). Instead of ensembles it 
is convenient to think of a (large) number of independent particles moving 
around in the sample. 



In a sample where the energy barriers are not identical, equation |5J 
does not hold at all time scales. Picture for example a 1-dimensional sample 
with mostly small energy barriers, E sma u, and a few much larger energy 
barriers, Ei arge , (see eg. [72|). At small time scales most of the particles will 



"think" they are living on an ordered sample with only the small energy 
barriers (since they haven't yet encountered a large energy barrier), and 
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the ensemble will follow equation 5.1 with a large diffusion coefficient (oc 



T(E sma ii)). However, at long time scales the particles will start to "feel" the 
effect of the large energy barriers, which will slow down the diffusion. At 
very long time and length scales the sample will appear homogeneous and 
the system again become diffusive, but now with a small diffusion coefficient 
(oc T(Ei arge )). 

The deviations from equation |5.1| can be quantified using the time de- 
pendent diffusion coefficient JT3|, D(t), or the frequency dependent diffusion 
coefficient [74|, D(s) (with s being the Laplace frequency, s = iu): 

D(t) = \j t {^X\t)) (5.2) 

poo 2 poo 

D{s) = s / e~ st D(t)dt = 1- / e- st (AX 2 (t))dt (5.3) 
Jo 2 Jo 

Note, that for diffusion in a sample with identical energy barriers, equation 



5J] leads to D(s) = D(t) = D. 

In a disordered sample one in general finds that the system becomes dif- 
fusive at long time scales (corresponding to low frequencies) when particles 
are moving much longer than the appropriate correlation length (assuming 
that such a correlation length exists); Dq = D(s — * 0) = D(t — > oo). In one 
dimension one finds D = a 2 (r -1 ) -1 , where the average is over the dis- 
tribution of jump rates, T. In higher dimensions no such simple expression 



exists, but for (3 — > oo one finds [ 76[ : -Do oc T(E C ), where E c is the so-called 
percolation energy (see section |5.3.2p . 

At infinitely short time particles never jumps more than once, and one 
finds (where Tl and Tr are respectively the jump rates to the left and to 
the right of a given site): 

{AX 2 (At)) = (a 2 r L At + a 2 T R At) (5.4) 
= 2a 2 (r)Ai (5.5) 

In the high-frequency limit, we thus find: = D(s — ► oo) = D(t —►()) = 
a 2 (T). 

In |2f| and paper I, the symmetric hopping model was treated as a 
model for frequency dependent conduction in glasses. In that context the 
particles represent non- interacting charge carriers (ions or electrons), and 
the lattice sites represents the positions in the glass where the charge carriers 
can reside. The frequency dependent conductivity, cr(s), is related to D(s) 



by the generalized Einstein relation [74|: 
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where n is the density of particles, and e is the charge of the particles. 

In the present work we shift the attention slightly, and treat the sym- 
metric hopping model more generally as a model for diffusion in disordered 
media. We focus on the "extreme disorder limit", i.e. low temperatures, 
where the model is found to exhibit universal behavior; D{s) becomes in- 
dependent (suitably scaled) of the temperature and the chosen probability 
distribution for the energy barriers, p{E). In the extreme disorder limit, it is 
not feasible to simulate the model using standard Monte Carlo (MC) tech- 
niques pi, 26, 77 1; Due to the large difference in jump rates (oc e~@ E ), the 
particle will jump many times over small energy barriers, and only "sample" 
the rest of the lattice on much larger time scales. Although this reflects the 
physics of the model in the extreme disorder limit, it makes MC methods 
unsuitable in this limit. 

In the following we set the attack frequency, To = 1, the lattice constant, 
a = 1, and Boltzmann's constant, ks = 1, thus defining the scales for time, 
length, and energy respectively. 



5.2 The Master Equation 

The starting point for both the analytical and the numerical methods, is the 
master equation jf78|-|82|] for the model; Let P(i, t\j, 0) denote the probability 
of the particle being at site i at time t, given that it was at site j at t = 0. 
The master equation for the system is then: 

dP{i ^ 0) = -liPii, t\j, 0)+J2 r(* - i)P(k, t\j, 0) (5.7) 

k 

where 7, = T(i — » k). The first term on the right hand side is the prob- 
ability flow out of site i, and the second term is the flow into site i. Defining 
the (N^xN^) matrix P(t) with the components Pij(t) = P(i,t\j, 0), we 
write the master equation on matrix form: 

|p(t) = HP(t) (5.8) 

H is here a (N T> xN' D ) matrix containing the jump rates; H# = — 7$ and 
Hj^fc = T(k — > i). Note that only TD + 1 elements in each row in H 
are different from zero (the diagonal and the TD elements corresponding to 
nearest neighbors); H is sparse. This is an essential feature in the numerical 
methods; without it we wouldn't be able to treat large enough sample sizes. 
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Taking the Laplace transform of equation 5.8 we get: 



sG(s)-P(t = 0) = HG(s) (5.9) 

where the Gjj(s) is the Green's function, i.e. the Laplace transform of 
Pjj(t): Gij(s) = J °° Pij(t)e~ st dt. The initial condition, P(t = 0), is given 
by the identity matrix, I, and we thus find: 



(sI-H)G(s) = I & G(s) = (si- H)- 1 (5.10) 

Before proceeding to discuss the analytical approximations and the nu- 
merical methods, we present a "preview" of the dynamics of the symmetric 



hopping model in figure 5.2. A 2 dimensional sample was set up using a 
box-distribution of barrier energies; p(E) = 1,0 < E < 1. The time evo- 
lution of P(i) was determined by a discrete version of the master equation 
(eq. P): 

P(t + At) = P(t) + AtHP(t) = (I + AtH)P(t) (5.11) 
A time step of At = 0.1 was used, and as initial condition a particle was 



placed at a particular site. In figure 5.2a P(i = 2) is shown for (3 = 0, i.e. 



infinitely high temperature. In this limit all the jump rates are identical 
(equal to unity), and as expected the probability distribution spreads in a 
symmetric manner. In figure |]|b P(t = 10 5 ) is shown for (3 = 40, with 
the same starting position, and same energy landscape (i.e. same set of 



energy barriers) as in figure 5.2a. The time at which the two temperatures 
are compared was chosen so that the probability of finding the particle 
at the starting point was the same (~ 0.04). The qualitative difference 
between the dynamics at high and low temperature is evident; At the low 



temperature (figure 5^2 b) the spread of the probability distribution P(i) is 
highly irregular, as a consequence of the particles "preferring" low energy 
barriers and avoiding high energy-barriers. The particular structure of P(t) 
at the low temperature depends strongly on where the particle was started. 
If it was started in one of the two enclosed empty sites (white in figure 
|5.2| b) it would be stuck there on the time scale used here, since these must 
be connected with high energy-barriers to the surroundings (otherwise they 
would not be empty in figure |T 
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Figure 5.2: P(i) for a 2 dimensional realization of the symmetric hopping 
model, with a box distribution of energy-barriers. Each site is represented 
by a square, and this is colored a shade of grey according to the value of 
P(t) at the given site, a) P(t = 2) for = 0. b) P(t = 10 5 ) for /3 = 40. 
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5.3 Analytical approximations 

In this section we'll briefly describe the 3 analytical approximations which 
will be compared to the numerical results, in section 15.61 and 



5.3.1 Effective Medium Approximation (EMA) 

In EMA [^,|8(| the disordered sample is replaced by an ordered sample (the 
"effective medium"), where all the jump rates are replaced by an "effective" 
jump-rate, Te(s). The value of Te(s) is determined by a self-consistency 
condition; A single jump rate in the effective medium is replaced by its value 
in the disordered sample, and an average is performed over the distribution 
of jump rates. The condition that this average should be equivalent with 
the effective medium leads to the following condition [25,|32], 80||: 



D-T \ 

=- ) = 5.12 

VD-{D-T){l-sG)/ T 

where G is the diagonal element of the Green's function for the effective 
medium, which depends on the spatial dimension, T>. In the extreme disor- 
der limit {(5 — > oo) one finds |^] for T> > 2: 

Dhx{D) = S (5.13) 

where D = D(s)/Dq, and s = s/s c , where s c is a suitably defined character- 
istic frequency. The EMA thus predicts universality in the extreme disorder 
limit; Properly scaled D{s) becomes independent of temperature and the 
distribution of energy barriers, p(E). 



In [25 1 the predictions of EMA was compared with numerical results 
using the ACMA method (described in section |5.4| ) in 2 dimensions. The 
existence of universality was confirmed by the numerical results (for 4 differ- 
ent energy-distributions), but the shape of D(s) was found to deviate from 
the one predicted by EMA. Or to be more specific: EMA predicts a shape 
of the universal D{s) that is growing to rapidly at the onset of frequency 
dependence; it is to "sharp" (figure 5 in f25|). The EMA prediction for 

-2 



{AX (t)) is discussed in 

5.3.2 Percolation Path Approximation (PPA) 

In the extreme disorder limit, the low-frequency diffusion is expected to be 



dominated by percolation |77, 82, 84,p5|; when diffusing over long distances 
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(i.e. in the low-frequency limit) the particles "chose" to do so by jumping 
over energy-barriers that are as small as possible. The largest energy barrier 
that a particle has to cross to move through the sample is given by the 



"percolation energy", E c , defined by |25|: 



p(E)dE=p c (5.14) 

o 

where p c is the percolation threshold for bond percolation (for an introduc- 
tion to percolation theory, see f8(|). In a 2 dimensional regular lattice p c 
equals 1/2 (exact) and in a 3 dimensional cubic sample p c equals 0.2488 (ap- 
proximated). Equation |5.14| can be interpreted as follows; In a given sample 
mark the energy-barriers (bonds) starting with the smallest energy-barrier, 
then the next smallest etc. At some point the marked energy-barriers will 
"percolate" i.e. make a cluster that stretches through the whole sample (the 
"percolation cluster"). The highest energy needed to get percolation is (for 
an infinite sample) E c . For the box distribution of energy barriers used in 
the present work (p(E) = 1, < E < 1), the percolation energy equals the 
(bond) percolation threshold: E c = p c . The percolation cluster is a fractal, 
with fractal dimension 1.9 and 2.5 in 2 and 3 dimensions respectively ]86| ]. 

In paper I we argue, that the reason why EMA does not predict the 
shape of the universal curve D{s) in the extreme disorder limit correctly 
might be, that it replaces the disordered sample by an effective media of 
the same dimension, whereas the actual low frequency diffusion happens 
on a cluster of lower dimensionality. We will here term this cluster "the 
diffusion cluster", and its fractal dimension Df (this is not the same as the 
percolation cluster, see below). PPA can be considered as a first attempt at 
trying to incorporate the fractal dimension of the diffusion cluster into an 
analytical approximation. 

The percolation cluster contain "dead-ends", which contributes little to 
the low frequency diffusion. Removing dead-ends from the percolation clus- 
ter leaves us with the "backbone", with fractal dimension 1.6 and 1.7 in 2 
and 3 dimensions respectively []86|]. In paper I we argue, that Df is even 
lower, since the backbone contains loops, where one of the branches usually 
will be preferred (since the jump-rates depend strongly on the energy barri- 
ers in the extreme disorder limit). In PPA the extreme view that Df = 1 is 
taken, i.e. that the (low frequency) diffusion takes place on 1-dimensional 
"percolation paths" . With this assumption we in paper I arrive at the PPA 
approximation, given by: 
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(5.15) 

PPA thus predicts universality in the extreme disorder limit, like is the 
case for EMA. In paper I , we find that PPA agrees better with the numerical 
results than EMA. This is especially true in 3 dimensions. The agreement is 
however not perfect, and a number of problems remain unanswered (see last 
section in paper I). The PPA prediction for (AX 2 (t)) is discussed in |87j. 

5.3.3 The Diffusion Cluster Approximation (DCA) 

The Diffusion Cluster Approximation can be thought of as a refinement of 
PPA (paper IV and V); instead of setting Df = 1, this is in DCA left as a 
parameter in the model. Using the approach of EMA, but with an fractal 
effective medium with fractal dimension Df, one finds: 

\Df/2 



\n{D) [Ij (5.10) 



This expression is limited to 1 < Df < 2. For Df > 2 DCA reduces to the 
EMA prediction (eq. 5.13| ), and for Df < 1 it is undefined. Since we do not 



have an independent estimate of Df we will in the present work treat it as 
a fitting parameter. This must obviously be taken into consideration, when 
comparing with EMA and PPA, since these have no fitting parameters. 

From the arguments given above, we expect Df to be limited above by 
the fractal dimension of the backbone. Furthermore we expect Df to limited 
below by the fractal dimension of the so-called "red bonds" , which are the 
singly connected bonds on the backbone, i.e. if a red bond is removed the 
backbone is broken into 2 parts. The fractal dimension of the red bonds are 
3/4 and 1.14 in 2 and 3 dimensions respectively. 
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5.4 The ACMA Method 

In this section we briefly describe the ACMA method, which was used to 
obtain the numerical results reported in [25| and paper I. 

Defining Aj as the x-coordinate of site i we can write the frequency 
dependent diffusion coefficient (equation |5.3|): 



„2 rco i 

D x (s) = - J e -*_^(^-A,) 2 (P(i,t|i,0))^ (5.17) 



2N d 



hi 



*i) 2 (G(i,*|j)) 



(5.18) 



'•J 



The subscript x in D x (s) is here used to emphasize that the diffusion coeffi- 
cient here is calculated from the mean square displacement in the x-direction. 
In more than 1 dimension similar expressions obviously hold for other direc- 
tions. 

Equation 5.10 and equation |5.18 together constitutes a method for com- 



puting D(s); calculate (G(i,s\j)} by inverting (si — H) (equation 5.10 ) and 
insert the result in equation |5.18| @]. 

r(i,i+D x 






N-l 



N 



Figure 5.3: The ac Miller-Abrahams (ACMA) electrical equivalent circuit 
for the symmetric hopping model. The admittance of the resistor between 
two sites is equal to the corresponding jump rate. The admittance of the 
capacitors are s. 



In my master thesis ]32|, |33| it was demonstrated, that the method de- 
scribed above is equivalent to the method of the ac Miller- Abrahams (ACMA) 
electrical equivalent circuit pBfl. In this method a large network of linear 



electric components are set up (see figure ll3 for an illustration of the 1- 
dimensional case). After eliminating all the 'internal sites' (those without 
numbers in figure |5.3| ) by the so-called generalized Star-Mesh transforma- 
tion, a(s) (oc D(s), see equation 5^) is calculated as a weighted sum over 
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the effective conductances between the 'external sites' (those with numbers 
in figure |5,3| ). The elimination process described above is equivalent with 
Gauss- Jordan elimination on the matrix (si — H) (3*|, and the Fogelholm 
algorithm [88] for the order in which the internal sites are eliminated is 
equivalent with the minimum degree pivoting algorithm JS9[j : At each step 
eliminate the site/row with the smallest number of connections/elements. 

A number of 'tricks' is used in the ACMA method, to improve efficiency 
over a straight forward computation of first equation |5.10| and then equa- 
tion |5,18| . The most important trick is that the (dense) matrix G(s) is 
not explicitly calculated, but instead a 'condensed' version where columns 
corresponding to the sites with the same x-coordinate are added together. 
Thus instead of solving for right-hand sides (inverting si — H) only N 
right-hand sides are used. 

In the simulations using the ACMA method presented in |25|] and paper 
I 'perfect electrodes' were used in the x-direction (see figure |5.3| ). This 
ensures that D(s) goes to a constant value as s — ► 0, corresponding to 
(AX 2 (t)) being proportional to t, as t — > oo, i.e. the system being diffusive. 
However, with this modification the model does not correspond directly to 
the solution of a master equation. This is evident from the fact that the 
mean square displacement as calculated in equation 5.17 is finite for any 
finite sample, and thus it can not be proportional to t at long times. 

To what extent the 'perfect electrodes' mimics the effect of real electrodes 
when experimentally measuring the frequency dependent conductivity, cr(s), 
will not be discussed here. What we are after here, is to calculate the bulk 
value of the frequency dependent diffusion coefficient -D(s) in the extreme 
disorder limit. 



5.5 The Velocity Auto Correlation (VAC) method 

In this section the Velocity Auto Correlation (VAC) method is developed. 
Its main advantages over the ACMA method is, that it can be (and is) 
used with periodic boundary conditions, and still give a diffusive regime 
for s — > (i.e. t — ► oo). Another advantage over the ACMA method 
(as it was implemented in the previous work), is that the VAC method is 
formulated in terms of sparse matrices, which means that standard methods 
for solving these can be usedR The differences between the two methods 



will be discussed further in section 5.5. 



lr The numerical results presented here was done using Matlab version 5.3.0. 
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To derive the VAC method we express the diffusion coefficient in terms 



of the velocity auto correlation function [74]: 



D x {s)= / (v(t)v(0)) x e- st dt (5.19) 
Jo 

v(t) is here the velocity (in the x-direction) of the particle at time t. The 
motion of the particle in the symmetric hopping model is instantaneous; 
when it jumps from one site to another, this happens in a infinitesimally 
small time-interval, At. We may thus assign the constant velocity ±a/At 
(where we briefly reintroduce the lattice constant, o) to the particle in the 
time interval, At, ensuring that it moves one lattice constant, a, either to 
the left or to the right. v(t) thus takes on the value a/ At in a time interval 
At when the particle jumps to the right, the value —a/At in a time interval 
At when the particle jumps to the left, and zero the rest of the time (It is 
not important whether or not this is a physically reasonable choice for v{t). 
The only requirement for eq. [5 . IS to hold is that x(t) = f o v(t)dt + x(o)). 
With this choice of v(t), the function v(t')v(0) has the value a 2 / At 2 if the 
particle jumps in the same direction at t = t 1 and t = 0, the value —a 2 / At 2 
if the particle jumps in opposite directions at t = t' and t = 0, and zero 
otherwise (i.e if the particle does not jump both at t = t' and t = 0). 

The probability, Prr, that the particle jumps to the right both at t = t' 
and t = is (for now we assume that t' ^ 0): 

hi 

At 2 

= xiT,TUj)(p(i,t\mm) (5.2i) 

hi 



Equation 5.20| can be read from the left as follows; If the particle starts at 



site j the probability of jumping to the right at t = is Tji(j)At. This 
means that it is now at site j ' + 1, and the probability that it moves to site % 
in the time t is {P(i, t\j + 1, 0)) , and finally the probability that it jumps to 
the right from site i is Tn(i)At. In equation p. 21 j is substituted with j — 1, 



and Tft(j — 1) = Ti(j) is used. Calculating in the same way the probability 
of the other events that contributes to the velocity auto correlation function, 
(v(t)v(0)) x (Prl, Plr, and Pll), we can now write (v(t)v(0)) x , in terms of 
the Green's function (where again use a = 1, and the delta function takes 
care of the special case t = 0, see below): 
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(v(t)v(0)) x = C5(t) + ^^2(P RR + P LL -P RL -P LR ) 



(5.22) 



(5.23) 



'■j 



[T L {j)T R {i) + T R (j)T L (i) - T L (j)T L (i) - T R (j)T R (i)} 
C5(t) + (5.24) 

l Td E( r «« " r ^*)) W, o)) (TlO") " r «(i)) 



Taking the Laplace transform of equation 5.24 we get: 



In the high- frequency limit, s — ► oo, we find from equation |5.10 
s|j)) — > 0. We thus get C = = (V) (see equation ^5| ): 



In the case where all the jump-rates are identical, T R (i) = Tl(i), equation 



5.26| leads to D x (s) = (V) = D^, as expected (see section |5.ip . 

Equation |5.26| can be used together with equation 5.10| to calculate 
D x (s). However this method poses a serious numerical problem at low tem- 
peratures where Do is very small compared to £>oo- Calculating D x (s) at low 
frequencies from equation 5.26| in this limit amounts to calculating a small 
difference between two large numbers, which leads to large uncertainties in 
the result. 

In the following we will derive a version of the VAC method, which does 
not have the numerical problems described above. We do this first in 1 



dimension (section |5.5.1 ) , followed by a general derivation in T> dimensions 
(section ^.5.2|) . 



5.5.1 The VAC method in 1 Dimension 



The idea behind the derivation of the VAC method in 1 dimension is to 
rewrite the problem in terms of the variable: 
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J R (i, t\j, 0) = r fl (i) [P(i, t\j, 0) - P(i + 1, t|j, 0)] (5.27) 

which can be interpreted as the particle current to the right of site i at time 
t given that the particle started at site j at t = 0. The master equation 
(|5.7|) can then be written as: 



dP( ^ |j ' 0) = J R (i ~ l,t\j, 0) - J R (i, t\j, 0) (5.28) 
at 

We now define two new matrices; T R is a diagonal matrix containing the 
jump rates to the right of each site: (T R )ij = 6(i,j)T R (i). A is a matrix 
with +1 on the diagonal, —1 on the sub-diagonal and zero everywhere else: 
(A)jj = 5(i,j) — 5(i — In more general terms; equals —1 if site j 

is the left neighbor to site i. Here (as in the rest of this chapter) periodic 
boundary conditions are implicit, i.e. (A)uv = — 1- We can now write the 
1 dimensional master equation as: 



-AJ fl (t), 3 R (t) = T R A T P(t) 



(5.29) 



These equations provides us with an an explicit expression for the matrix 
H in the master equation ( |5.8j ): H = —AT R A T . What we are after here is 
however an equation for J#. To get this we take the Laplace transform of 
equation 5.2£: 



sG(s) + AJ R (s) =P(i = 0) 



(5.30) 



where Jr(s) is the Laplace transform of Jfi(t). By multiplying from the left 
by A T and using Jr(s) = TrA t G(s) we get: 



A T & 



(sT^ 1 + A T A) 3 R (s) 
3r(s) = ( S T~ 1 + A t A)- 1 A t 



(5.31) 
(5.32) 



Here we assume that the diagonal matrix Tr is invertible, i.e. that all 
the jump rates are different from zero. We are now ready to derive an 
equation for D(s) in 1 dimension. To do this we rewrite equation 5.26 using 
T L (i) = T R (i - 1): 



D(s) = (T R ) 



*,3 



F R (i - 1)) <G(i, s\j)) (T R (j) - F R (j - 1)) 
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This can be written as (where 1 is a column vector containing all l's): 



D(s)N 



i t t r i - i T r R A T G(s)Ar R i 



l 1 T R 1 - l 1 J R (s)AT R l 
1 T T R 1 - l T (sT^ + A T A)" 1 A T Ar i jl 
^{sT^ + A T A)- 1 [(sT' 1 + A T A) - A T A] TrI 
sl^sT" 1 + A T A)- 1 1 
sl T x, (sT~ l + A T A)x = 1 



(5.33) 
(5.34) 
(5.35) 

(5.36) 
(5.37) 



We have now avoided the problematic subtraction (in equation 5.26| and 
5.33| ), and reduced the problem to finding the vector x by solving a sparse 
system of linear equations. 



5.5.2 The VAC method in V Dimensions 

In this section we generalize the VAC method to T> dimensions. In analogy 
with JR(t) used in the 1-dimensional version, we define Jfe(t) as the particle 
current in the fe'th direction: 

J fc (i) = r fc H£P(t) (5.38) 

where Tk and are generalization of the matrices Tr and A used in the 
previous section. r& is a diagonal matrix, which for each site contains the 
jump rate to the "right" in the fe'th direction. (H^)^- equals —1 if site j is the 
"left" neighbor to site i along direction k, and it has 1 on the diagonal and 
zero everywhere else. The structure of depends on the numbering scheme 
chosen for the sites, but the resulting physics is obviously independent of 
this. An explicit expression for (and the corresponding numbering of 
the sites) is given in appendix A. 

The goal is now to set up a master-equation for the currents, and from 
that calculate D(s), like it was done in the 1-dimensional case. Obviously 
it is not possible to solve for Jfc(s) without at the same time solving for the 



currents in the other directions. We thus generalize equation |5.38| to 

J*(t) = r*H^P(t) J*(s) = r*H^G(s) (5.39) 
where we have defined the following block-matrices: 
/ Ji(t) \ / Ti \ 



J*(t) 



J 2 (t) 



o r 5 



o o 



o 



H* = (Hi, . . . , H 



(5.40) 
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and J*(s) is the Laplace transform of J*(i). 

The master equation for P(t) is similar to the one in 1-dimension (eq. 
5.29| ) , but the change in the probability at a given site now has contributions 
from all T> directions: 



d 



v 



-P(t) = - H fc J fc W = -H*J*(i) = -H*r*H* P(t) O (5.41) 

fc=i 

sG(s) +H*J*(s) = P(t = 0) (5.42) 

By multiplying equation |5.42| from the left with r*H^, and substituting 
for J*(s) (equation 5.39 ), we arrive at the (Laplace transformed) master 
equation for J*(s): 

(sI + I%HrH*) J # ( s ) =r,HrP(* = 0) & (5.43) 

J*(s) = (si + r.H^H*) -1 r*H* P(* = 0) (5.44) 

To derive the final result for Dk(s) (i.e. the diffusion coefficient as calcu- 
lated from velocity correlations in direction k) we define as a column vec- 
tor with T>N V elements, where the N v elements corresponding to currents 
in the fc'th direction is 1 and the rest is zero (so that eg. l^r^l^ 

/0\ 



(5.45) 



V / 



From equation 5.26| we finally get (using P(i = 0) =1): 



D k (s)N 



v 



r fc - j fc (s)H fc r fc i 



i— j*(s)h* r*i 



* ■ L fe 



(5.46) 
(5.47) 
(5.48) 



\\ (\ - (si + r*H^H*) 1 r*H^H^ r*i k 
i% (si + r^H^H*) ^ ( (si + r^H^H*) - r*H^H*^ r*i* 



sll" (si + r*H^H*) ^,1* 

si^x , (sr; 1 + h^hJ x 



(5.49) 
(5.50) 
(5.51) 
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Like in the 1-dimensional case in the previous section, we have now re- 
duced the problem to finding the vector x by solving a sparse system of linear 
equations. In T> dimensions the matrix involved in this is a (VN T 'xT)N v ) 
matrix. 

For real s > the matrix (sT" 1 + H^H*) is positive definite^], which 
ensures the numerical stability when doing Gaussian elimination [p9j]: 

y T (sT: 1 + H*) y = sy^y + y T H^H*y (5.52) 

= sy T r- 1 y + (H*y) T (H,y) (5.53) 

The second term is greater than or equal to zero, and the first term is 
greater than zero, since T" 1 is a diagonal matrix with only positive values 
on the diagonal. When representing the problem in finite precision in the 
computer this last statement is violated for very low s values; The elements 
corresponding to small energy barriers (large jump-rates) are effectively set 
to zero when the two terms are added (all elements in H* are 0(1)). This 
means that the matrix in practice is not positive definite and the attempt 
to solve it fails. To avoid this numerical problem a small constant 5 (= 
10~ 14 ) is added to the diagonal ensuring that the matrix is positive definite 
and thus can be solved. In physical terms this corresponds to applying 
a minimum value for the energy barriers, ensuring that all jump-rates are 
finite (< s/S). It was checked numerically (by varying 5) that the effect of 
this procedure is negligible (relative error < 10~ 4 ), as expected from our 



physical understanding of the model (see section |5.3.2|) . 

Like in paper I the jump rates corresponding to energy barriers larger 
than E c + K/j3 (K = 6.4) was set to zero, to speed up the calculations (and 
reduce the amount of memory needed). By varying K the relative error 
introduced by this procedure was estimated to be less than 1%. 



Equation |5.51 was solved using Cholesky factorization for real s-values 
and LU factorization for imaginary s- values. In both cases pivoting was 
done using the minimum degree algorithm, and all computations was done 
using Matlab version 5.3.0. 



2 The matrix A is symmetric and positive definite if y T Ay > for any y. 
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5.6 Numerical results, D(s) 



In this section we report the results of numeric calculations of D(s) in 3 
dimensions using the VAC method. Like in and paper I the calcula- 
tions are done with for real Laplace frequencies, s, i.e. corresponding to 
imaginary values of ui [s = ioS). When comparing the numeric results and 
analytical approximations for real values of s, we are relying on the principle 
of analytical continuation |kJ; If two complex functions coincide on a line 
in the complex plane, they coincide everywhere in the complex plane (where 
they are both well-defined). 



In section |5.7| we present the corresponding results for real values of u>. 
Besides providing results directly comparable with experimental results for 
cr(uj), it will be clearly demonstrated that we can indeed trust the principle 
of analytical continuation when comparing analytical approximations and 
numerical results. 



It has already been demonstrated in [|25j] and paper I that the symmetric 
hopping model becomes universal in the extreme disorder limit, i.e. that 
-D(s) (or equivalently 5"(s)) becomes independent of temperature and the 
distribution of energy-barriers. Here we will focus on the shape of D(s), 
and only present results for the box-distribution of energy barriers: p(E) = 
1,0 < E < 1. 



Figure 5.4 show the frequency dependent diffusion coefficient for real 



Laplace frequencies, D(s), in a log-log plot, for 4 /3-values. Both axis in 



figure 5.4 are scaled by the DC level, Dq = D(s — > 0). In the inset of figure 



5.4 is shown Dq scaled by T(E C ) vs. f3. For > 80 Dq is seen to be well 



approximated by: 



D cx /T 7 r(E c ), 7 = 0.89 ±0.01 (5.54) 



Note that the dominant /^-dependence of Dq is given by T(E C ), which changes 
by 30 orders of magnitude for the /3-values used (see table |5.1| ) . In the sim- 
ulations reported in paper I , we found: 7 = 0.81 db 0.04 |32f. 
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1.2 • 10~ 1U 


1.5 • 10~ iy 


4.2 • 10~ 3Y 


Sc 


6.3 • l(T uti 


7.1 • 10" 11 


3.4 • 10~ 2U 


3.7 • 10-' M 



Table 5.1: Various parameters as a function of the /3-values used. The results 
reported are from averages over 100 NxNxN samples, where N depends on 
(3 as shown above. 




-6 -4 -2 2 4 6 8 10 



log 10 (s/D ) 

Figure 5.4: The frequency dependent diffusion coefficient, D(s), vs. the 
(real) Laplace frequency, s, both scaled by Dq. Inset: Dq/T(E c ) vs. (3, a 
power-law fit, and a power-law with exponent —1 for reference. 

To illustrate how D(s) approaches a universal curve (suitable scaled) at 
high /3-values, the data presented in figure ^4] was scaled by the charac- 
teristic frequency, s c , defined here by D{s c )/Dq = ylO. In figure |5.5| the 
result of scaling the data in this way is shown. It is clearly seen that D{s) 
approaches a universal curve, as the /3-values increases. Or in other words: 
the way the system approaches (long time) diffusion becomes universal as 
the temperature is decreased. The universal curve is estimated to be close 
to the data for (3 = 320, and the frequency regime for which the data follows 
the universal curve is seen to increase as j3 increases. 
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In the inset of figure [53] is shown s c scaled by Dq vs. (3. For (3 > 80 s c 
is seen to be well approximated by: 



Scfxp-^Do, 7 = 1.37 ±0.01 



(5.55) 



The dominant /3-dependence of s c is given by Dq, which changes by 31 
orders of magnitude for the /3-values used (see table |5.1| ) . In the simulations 
reported in paper I, the power law was less well defined, and rj was estimated 



to be 1.48 db 0.06 (depending on the /3-range used in the power-law fit) [32]. 



i2- 



o 




2 

log 10 (s/s c ) 



10 



Figure 5.5: Same data as in fig. but scaled on the frequency axis, to agree 
at log 10 (D(s c )/Dq) = 0.5 Approach to universal curve is evident. Inset: The 
scaling parameter, s c divided by Dq vs. (3, and a fit to a power-law. 



For both of the scaling parameters, Dq and s c , the /3-dependence given 
by the power-laws reported are very small compared to the /3-dependence 
given by the other involved quantities. In the simulations presented here, 
the main focus has not been on determining the scaling exponents, 7 and r/, 



precisely, but on the shape of the universal curve seen in figure 5.5 



Before proceeding to compare the numeric results with the analytical ap- 
proximations, we will briefly discuss the possibilities of computing (AX 2 (t)) 
from D(s). As discussed above we will compare numerical results and analyt- 
ical approximations for D(s), and consequently we are calculating (AX 2 (t)) 
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Figure 5.6: Estimating (AX 2 (t)} from D(s). a) Fit to D(s) to functional 
form so that 2D(s)/s 2 (inset) has known inverse Laplace transform b) Re- 
sulting estimate for (AX 2 (t)), and (in inset) D(t) and (AX 2 (t))/2t (g(x) is 
here the gamma-function, normally denoted T{x)). 
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only for illustrative purposes. From equation 5.3 this is in principle straight 
forward; (AX 2 (t)} is the inverse Laplace transform of 2D(s) /s 2 . The general 
inverse Laplace transform is however in practice problematic |9(|, and here 
we furthermore have the problem of the many decades of frequency /time 
involved. Instead we use the following trick; we fit D(s) in such a way 
that 2D(s)/s 2 has a known inverse Laplace transform. This is done for 
(5 = 320 in figure |5.6| a, where D(s) is fitted to a sum of power-laws. Note 
that this is a purely empirical fit, and it does not necessarily have anything 
to do with the true functional form of D(s). The inset in figure [5lj| a shows 
2D(s)/s 2 and the corresponding fit, i.e. the quantity to which the inverse 
Laplace transform is applied. The resulting (AX 2 (t)) is shown in figure 
|5,6jb, without using the scaling parameters, so that the actual size of the 
quantities is seen (t c = l/s c ). Note that the system becomes diffusive at 
very long time scales; t « 10 40 . This illustrates the impossibility of using 
traditional MC techniques, since this would require at least 10 40 time steps. 
The inset in figure 5.6b shows D(t) and (AX 2 (t)/2t) as resulting from the 
procedure described above. 

In figure |T7] the universal curve for D(s) (represented by the data for j3 = 
320) is compared with the 3 analytical approximations described in section 
|5.3| . The fractal dimension Df used in the Diffusion Cluster Approximation 
(DCA) was treated as a fitting parameter. The fitting procedure will be 



explained when discussing figure |5.9| . The fit of DCA is seen to be almost 
perfect, and clearly better than both PPA and EMA. This is perhaps not 
surprising given that DCA has a fitting parameter which PPA and EMA has 
not, and the value of fit achieved by DCA thus greatly depends on whether 
an independent argument can be given for the value of Df. For now we only 
note, that the value Df = 1.35 is in the expected range: 1.14 < Dt < 1.7 



(see section 5.3) 



Like in paper I, we now proceed to focus on the shape of D(s), by 
calculating the apparent power-law exponent, i.e. the logarithmic derivative 
of D(s): 

ologjoW 

When plotting \i as a function of log 10 (-D(s)/.Do) we can compare the shapes 
of the data presented above, without relying on an empirical scaling of the 



frequency axis. This is done in figure |5,8j , for the same data as shown in 
|5.5| . Plotting the data in this way is seen to be more sensitive; The small 
difference between the data for j3 = 160, and j3 = 320 in figure |5T^ is more 
pronounced in figure 
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Figure 5.7: Comparing the universal curve for D(s) {(3 = 320), with the 3 
analytical approximations. The approximations were scaled to fit the data 
at the highest frequencies. The fractal dimension, Df = 1.35 in the Diffusion 
Cluster Approximation (DCA) was determined by fitting (see fig. |5.9| ). 
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log 10 ( D(s)/D ) 

Figure 5.9: Comparing for the universal curve (approximated by (5 = 320) 
with the three analytical approximations. The fractal dimension, Df = 1.35, 
in DCA was fitted (by hand) to the numerical data in this plot. 



In figure |5.9| we compare the universal curve for [i (approximated by the 
data for j3 = 320) with the 3 analytical approximations. As in figure |5.7| the 
data from the simulations is seen to lie between EMA and PPA, and DCA 
is seen to give a good fit. The fitting of DCA was done by hand using this 
figure. The inset in figure |5.9| shows the same data, but focusing on the low 
values of D, corresponding to low values of § and /x. In this limit, the data 
is seen to deviate from DCA, and seems to approach EMA. 

The behavior of D(s) at low frequencies, i.e. the departure from the DC 
level (Do) is better investigated by plotting D(s) — Dq as its done in figure 



5.10| . The universality for the numerical data is seen to hold even at the 



very low frequencies in this plot. This demonstrates that the universality 
seen at low frequencies in D(s) (figure |5.5|) , is not just a consequence of 
the DC level being dominant. On the other hand the apparent reasonable 



agreement between D(s) and DCA seen in figure 5/7 breaks down when the 
DC level is subtracted. 

The low-frequency expansion of D(s) is known to be p2|| : 

D(s) = 1 + As + Bs 3/2 + ... (5.57) 



In figure 5.1C the numerical data and EMA is seen to agree with the first 



two terms of this expansion. 
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Figure 5.10: Focusing on the very low frequencies: D(s)/Dq— 1 vs. s/s c . The 
universality for the numerical data holds even at the very low frequencies. 
The numerical data and EMA are well approximated by D(s) = Dq(1 + As), 
thus agreeing with the first two terms of eq. 
PPA and DCA. 



5.57. This does not hold for 



5.7 Numerical Results, D{uS) 

In this section we present numerical results for D(u), i.e. for imaginary 
Laplace- frequency, s = iw. In this case the diffusion coefficient is a complex 
quantity, and we can compare both the real and imaginary part with the 
analytical approximations. 



In figure 5.11 is shown the (scaled) real part of D (u>) vs. the frequency 



(scaled). The scaling parameters, Dq and lo c is shown in figure ^3.13 , Like 



we found in the previous section for real Laplace frequencies the data is 
here seen to approach a universal curve which agrees well with DCA, with 
Df = 1.35. This is the value found in the previous section from D(s); it was 
not found to be necessary to make a new fit for D{uj). 

Figure 5.12| shows the imaginary part of D(oj) corresponding to the data 



in figure 5.11, and using the same scaling parameters (and Df). Universality 
is evident at low frequencies (agreeing with EMA, apart from scaling), and 
is approached at higher frequencies. 
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Figure 5.11: Real part of D(u) vs. u, both scaled (see fig. 5.13 ). Df = 1.35 
was used in DCA, like for real Laplace frequencies (fig. 5.9), i.e. it was not 
fitted to these data. 
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Figure 5.13: Scaling parameters used for data with imaginary Laplace fre- 
quencies, s = iu. Note the agreement with the scaling parameters used for 
real s (fig. jA, and fig. 5.5). 



As an analogue to the (apparent) exponent fj. (eq. |5.56|) used for D(s) in 
the previous section, we can define an exponent for the real part of D(u): 



l^real 



rilogioQP'M) 



(5.58) 



H rea i is plotted in figure 5.1^ as a function of cu. Note that the convergence 
toward universality is more "abrupt" than it was found for D(s) (fig. |5.8j ); 
Only the data for j3 = 40 deviates significantly from DCA. 



In figure 5.15 we focus on the low frequency behavior of D'{uj) by sub- 
tracting the DC level, which gives us the chance to check the agreement with 
third term in the low-frequency expansion (Eq. |5.57| ), i.e. the exponent 3/2. 
Neither of the approximations agrees with this, which in particular means 
that EMA only agrees with the first 2 terms in the low-frequency expansion. 
The numerical data seems to agree better with the exponent 3/2, although 
it is difficult to judge the numerical significance of this. 
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Figure 5.15: Focusing on low frequencies: D'(w)/Dq — 1 vs uj/uj c . Dashed 
thin curve: low frequency expansion (exponent 3/2, eq. 5.57 ). 
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We define the apparent exponent for the imaginary part of D(u) in the 
following way: 

Himaq = 7~. 7 (5.59) 

alogioM 



This is plotted as a function of D"{oo) in figure 5.1£. EMA has the right 



value at the very low frequencies, as seen in fig. 5.12. DCA has the right 



behavior from D"(uj) > Dq, with small deviations at the highest values 
(notice the y-axis being different from the other exponent-plots). 
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Figure 5.16: The exponent [limag (eq. |5,5S| ). 

In the data presented so far we have found much the same behavior in 
D(u) as for D(s) with regards to how the numerical data agrees with the 
analytical approximations; at high frequencies the numerical data falls be- 
tween EMA and PPA and is well fitted by DCA with Df = 1.35, whereas the 



very low frequencies is governed by the low frequency expansion (eq. 5.57 ), 
which (for the first 2 terms) agrees with EMA. In contrast the approach to 
universality seems to different for especially D(s) and D'{oj). The quantity 
that approaches universality in a similar manner as D{s) is the absolute 



value of D(u>). This is illustrated in figure 5.17 where we plot the apparent 
exponent for the absolute value of D(uj): 

dlog 10 (|DM|) 



H-abs 



d\og 



10 



(5.60) 
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Figure 5.17: The exponent n a b s (eq. |5.60j ). The approach to universality is 
almost identical to the one seen for real Laplace frequencies (fig. |5.9|). 



5.8 VAC vs. ACMA 

In figure |5.18 the universal curve for D(s) as calculated from the VAC 



method (fig. 5.7 ) and the ACMA method (fig. 1 in paper I) is compared. 
There is a small but significant difference between the results from the two 
methods, as can also be seen by comparing the apparent exponents from 



the two methods (fig. 5.8 and fig. 2 in paper I). The main difference be- 
tween the two methods is the boundary conditions, so the differences in 
these are the "main suspect" for the (slightly) different results. At infinitely 
large samples we would expect the results of both methods to converge to 
the bulk-limit. The way the two methods converge to the bulk limit might 
however be different; In the ACMA method there is some fraction of the 
sites that are close to an electrode, and thus give a wrong contribution to 
D(s). As the sample size increases this fraction decreases, and the results 
converges to the true bulk- limit. In the VAC method there are no sites that 
are "worse" than the other, and the finite size effects are more subtle; they 
arise from particles traveling trough the sample to where they started, and 
thus experiencing false correlations. 
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Figure 5.18: Comparing D(s) for (3 = 320 calculated by the VAC method 
(fig. |5.7| ) and the ACMA method (fig. 1 in paper I). A small but significant 
difference is found. 



In figure 5.19 we compare the apparent exponent // for f3 = 320 as 
calculated with N = 96 (i.e. the data in fig. [5.8D and N = 64 (averaged 
over 600 samples). The agreement is seen to be excellent, with the largest 
error at the very low frequencies. The inset in figure 5.19 shows the In- 
dependence of Do, with error-bars estimated from the standard deviation. 
The tendency to converge to a constant as N increases is evident. 



In figure |5.19 is shown /i calculated independently for each of the 100 
(96x96x96) samples at (3 = 320, i.e. without averaging. The noise shows 
that the samples are not self-averaging, i.e. it is necesarry to average over 
a number of samples. The data points are distributed evenly around DCA, 
showing that the averaging over (non self-averaging) samples do not intro- 
duce systematic errors (compare fig. 
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5.9 Conclusions 

The main result in this chapter is the development of the Velocity Auto 
Correlation (VAC) method. At first it might seem strange to work in terms 
of velocities in a hopping model, as it is evident from the following quote 
from Scher & Lax ||74| : 



The relation as it stands [ i.e. eq. [5.19 in the present work 



relating D(s) to the velocity auto correlation function ] is in- 
convenient to use in a hopping model since it refers to velocities 
rather than positions. 

It should be evident by now, that it is worth to suffer the (initial) incon- 
venience of working with velocities; The VAC method is clearly to prefer to 
the ACMA method, since it can be used with periodic boundary conditions, 
and still give a diffusive regime. The physical reason for this being possible 
is that the diffusive regime is characterised by loss of correlations, which 
is possible in a finite sample (as opposed to the mean square displacement 
being proportional to time). 

The numerical results was found to be in excellent agreement with the 
Diffusion Cluster Approximation, except for the very low frequencies where 
the low frequency expansion holds. The agreement with DCA was achieved 
by a (single) fit to the fractal dimension of the diffusion cluster, Df = 1.35, 
which is within the limits expected from percolation arguments: 1.14 < 
Df < 1.7 (see section |5.3| ). Two questions needs to be answered, to decide 
whether the agreement between the numerical data and DCA signals that 
the picture behind the approximation is correct, or if it is simply a result of 
fitting: 

• Is DCA a good approximation of diffusion on a diffusion cluster with 
fractal dimension Df? This corresponds to the check of PPA in 1 
dimension shown in figure 2a in paper I. 

• Is the diffusion cluster (in 3 dimensions) a fractal with fractal dimen- 
sion (close to) 1.35? 

If either of the answers to these questions are negative, then the agreement 
found between numerical data and DCA is a coincidence, and DCA is merely 
a convenient way to describe the data. It is left as future work to answer 
these questions. 
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Appendix A 

Numbering scheme for VAC 
method 



The sites in the P-dimensional regular lattice are numbered by: 



v 

SiteIndex = J](C i -l)iV^ 1 + l (A.l) 
3=1 

where Cj is the coordinate (counted from 1 to N) of the site in the j'th 
direction. 

The (N v x N v ) matrix H fc used in section |5.5.2| (eq. |5.38|) describes the 



connectiviy in the lattice; (H.k)ij equals —1 if site j is the "left" neighbor to 
site i along direction k, and it has 1 on the diagonal and zero everywhere 
else. With the numbering scheme given above, H/% is given by 1(8)1(8). ..A...® 
I (8 I, where I is the (N x N) identity matrix, (8 is the direct (Kronecker) 
multiplication, and the matrix A is at the j'th position from the right]]. 



1 The (JV x N) matrix A is here defined as in section 5.5.1: (A)y = 5(i,j) — S(i — 1, j) 
(with periodic boundary conditions) 
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APPENDIX A. NUMBERING SCHEME FOR VAC METHOD 



In 2 dimensions with N=3 we eg. have: 
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